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ABSTRACT 


The  buckling  analysis  of  imperfect,  thin,  circular,  cylindrical, 
stiffened  shell  under  uniform  axial  compression,  for  various  boundary  con- 
ditions is  first  investigated.  A methodology  is  presented  for  predicting 
critical  conditions  for  such  configurations.  This  methodology  is  based  on 
the  smeared  technique  and  the  von  KArmAn -Donne  11  nonlinear  kinematic  relations 
in  t lie  presence  of  geometric  imperfections.  The  computational  procedure 
employs  a Fourier  series  type  of  separated  solution  and  through  the  Galerkin 
procedure  the  field  equations  are  reduced  to  a system  of  ordinary  differential 
equations.  These  equations  are  solved  by  the  finite  difference  scheme. 
Numerical  results  for  numerous  stiffened  and  unstiffened  configurations  are 
presented . 

Then  the  effect  of  initial  geometric  imperfections  on  the  optimal 
stiffened  circular  cylindrical  shell  sunder  uniform  axial  compression  is 
assessed.  The  imperfection  sensitivity  of  geometries  corresponding  to  values 
of  the  design  variables  surrounding  the  optimal  configuration  is  investigated 
for  two  design  configurations.  A design  methodology  is  proposed  through 
which  one  may  arrive  at  the  minimum  weight  configuration  in  the  presence  of 
geometric  imperfections  of  predetermined  maximum  amplitude  and  of  a shape 
that  yields  the  greatest  reduction  in  the  linear  theory  buckling  load. 
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INTRODUCTION 

Stability  of  thin  circular  cylindrical  shells  (with  or  without 
stiffening),  because  of  its  importance,  has  enjoyed  tremendous  attention 
for  the  past  seventy  years.  Although  a complete  understanding  of  all 
the  details  of  the  phenomena  involved  has  not  yet  been  reached,  it  has 
been  well  established  that  the  discrepancy  between  classical  theoretical 
predictions  and  experimental  results  lies  primarily  in  the  fact  that  the 
system  is  sensitive  to  geometric  imperfections,  the  presence  of  which 
is  unavoidable. 

The  imperfection  sensitivity  of  the  system  was  initially  established 

through  strict  postbuckling  analyses  of  the  perfect  geometry  system.  In 

addition,  it  was  explained  that,  the  load  carrying  capacity  of  such 

systems  is  directly  related  to  the  lowest  load  corresponding  to  post- 

buckling  states  of  equilibrium.  The  first  theoretical  investigation  of 

this  type  was  reported  by  von  Karman  and  Tsien^  in  1941.  The  investigators 

calculated  postbuckling  equilibrium  states,  for  an  unstiffened,  axially 

compressed,  thin,  circular  cylindrical  shell,  corresponding  to  loads  far 

below  the  classical  critical  load.  The  calculations  were  based  on  a 

number  of  simplifying  assumptions  the  most  important  being  the  neglect  of 

2-4 

the  effect  of  the  boundary  conditions.  Many  subsequent  investigators 
attempted  to  improve  the  calculations  of  von  Karman  and  Tsien  in  order  to 
find  the  smallest  postbuckling  equilibrium  load.  This  search  came  to  an 
end  when  Hoff,  Madsen  and  Mayers'"  found  in  their  calculations  that  the 
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minimal  pos tbuck 1 ing  load  tended  towards  zero  with  improved  functional 
representation  (taking  more  and  more  terms  in  the  series)  for  the  solu- 
tion to  the  governing  equations  and  with  diminishing  thickness.  In 
addition,  Madsen  and  Hoff^  repeated  the  investigation  by  employing  more 
accurate  kinematic  relations  than  those  of  Donnell.  The  difference 

between  these  results  and  the  previous  ones  was  insignificant.  Futher- 
7 / / 

more,  Koiter  has  shown  that  the  von  Karman-Donnell  equations  are  also 
applicable  to  problems  with  arbitrarily  large  displacements,  when  the 
function  describing  the  radial  displacement  is  taken  to  be  the  curvature 
function  defined  in  his  paper. 

g 

Koiter  was  the  first  to  question  the  use  of  the  minimal  postbuckling 
equilibrium  load  as  a measure  of  the  load  carrying  capacity  of  the  configu- 
ration. Instead,  he  proposes  to  find  the  critical  load  (limit  point)  of 
the  imperfect  system.  Koiter's  work  is  the  first  attempt  to  the  buckling 
analysis  of  an  imperfect  shell,  but  his  method  is  limited  to  the  neighbor- 
hood of  the  classical  load  and  therefore  to  small  imperfections  of  certain 

spatial  form.  This  approach  has  been  adopted  by  many  investigators  in- 

9 

eluding  Hutchinson  and  Amazigo  who  treated  the  stringer  or  ring  stiffened 
thin  cylindrical  shell.  Excellent  reviews  on  the  subject  may  be  found  in 
the  works  of  Hoff^,  and  Hutchinson  and  Koiter^. 

Many  of  the  postbuckling  analyses  that  are  based  on  Koiter's  proposi- 
tion (see  Ref.  11)  disregard  the  effect  of  end  conditions  by  assuming  that 
the  cylinder  length  is  extremely  large.  A systematic  experimental  investi- 
gation dealing  with  cylinders  of  various  lengths  (Ref.  12)  revealed  that 


Lite  postbuckling  behavior  is  strongly  influenced  by  the  cylinder  length. 

1 3 

Narasinham  and  Hoff  ’ analyzed  an  unstiffened,  thin,  circular, 
cylindrical,  imperfect  shell  of  finite  length  under  uniform  axial  com- 
pression. They  solve  the  nonlinear  equations  by  employing  a separated 
series  solution  for  the  dependent  variables,  each  term  of  which  contain^', 

- P 

a function  of  x multiplied  by  a cosine  term  in  y (Fourier).  Thus  the 
equations  are  reduced  to  a system  of  ordinary  differential  equations, 
which  in  turn  are  solved  by  the  finite  difference  scheme.  Although  the 
equations  are  developed  for  arbitrary  terms  of  Fourier  series,  the  solu- 
tion is  restricted  to  just  one  term  for  the  displacement  function.  A 
similar  procedure,  but  one  that  employs  the  "shooting  method"  (Ref.  15) 
instead  of  the  finite  difference  technique,  is  employed  by  Arbocz  and 
Sechler^  In  their  investigation  of  the  buckling  behavior  of  axially 
compressed  imperfect  cylindrical  shells. 

With  the  exception  of  the  work  of  Ref.  9,  there  is  virtually  no 
reported  investigation  on  the  buckling  behavior  of  imperfect  stiffened 
configurations.  The  first  part  of  the  report  presents  a methodology  for 
analyzing  the  buckling  of  a uniformly  compressed,  stiffened  (rings  and 
stringers),  thin,  circular,  cylindrical  imperfect  shells  of  finite  length 
and  various  boundary  conditions.  The  analysis  employs  the  von  Karman- 
Donncll  large  displacement  equations  and  the  smeared  technique.  The 
solution  procedure  is  similar  to  that  of  Ref.  13  but  the  limitations  on 
the  spatial  character  of  the  imperfection  has  been  relaxed  considerably. 
Results  have  been  produced  for  special  case  geometries  that  have  been 
reported  in  the  open  literature  (bench  marks)  and  for  new  configurations 


r 


of  the  stiffened  type  in  both  directions. 

The  latter  part  of  the  present  report  applies  this  method  to  the 
investigation  of  the  effect  of  initial  geometric  imperfections  on  the 
optimal  (linear  theory)  stiffened  cylinder  configuration  under  uniform 
axial  compression. 


4 


CHAPTER  II 


MATHEMATICAL  FORMULATION  OF  THE  NONLINEAR  BUCKLING  ANALYSIS 

, By  employing  the  von  Karman-Donnell  kinematic  equations  for 

geometrically  imperfect  [w°(x,y)]  thin  cylindrical  shells  one  can  easily 
derive  the  compatibility  and  transverse  equilibrium  equations  in  terms 
of  the  radial  displacement  w and  the  stress  function  F,  as  well  as  the 
expressions  for  the  total  potential  and  the  "unit  end  shortening".  The 
procedure  employed  is  similar  to  that  outlined  in  Ref.  13  for  unstiffened, 
imperfect  thin  cylindrical  shells. 

Consider  a geometrically  imperfect  stiffened  cylinder  under  uniform 

o 

axial  compression.  Let  w (x,y)  denote  the  deviation  of  the  shell  mid- 
surface  (taken  to  be  the  reference  surface)  from  the  corresponding 
perfectly  cylindrical  one.  Let  u,v,  and  w denote  the  displacements  of 
material  points  on  the  reference  surface  (see  Fig.  1). 

The  kinematic  relations,  first  proposed  by  Donnell^  are  given  below 


e = u,  + ^(w^  + 2w,  w°  ) 

xx  x x x x 


e = v,  - w/R  + ^(w^  + 2w,  w°  ) 

yy  y y y y 


v =2e  = u,  + v,  +w,w,  + w,  w?  + w,  w° 

'xy  xy  y x x y y x x y 


(1) 


H =w>  ; K =w,  ; H = w, 

xx  xx  yy  yy  xy  xy 


The  stress  and  moment  resultants  to  strains  and  changes  in  curvature 
and  torsion  are  taken  from  Ref.  18.  They  were  derived  by  employing  the 
smeared  technique. 
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N - E [(1+X  )e  +ve  - e X h ] 

xx  xx  xx  xx  yy  x xx  xx 

P 

N = E [ve  + (1+X  )e  - e X h 1 

yy  xxp  xx  yy  yy  y yy  yy 


N - E [ 1-v) e 1 
xy  xx  xy 

P 


12  2, 

i 

t 


12 


M = D i (1+p  ) + -%  e“X  |h  + wt  - -w  e \ e f 

xx  ..  Kxx  x xxJ  xx  yy  2 x xx  xxJ 


(2) 


12  2,  "j  12 

) ' “T>  e A H ” ~ c /v  o | 

yy  t2  y yyJ  yy  fc2  y yy  yyJ 


M = Di  vx  + ( 1+p  ) + -?;  e A I k - e A e 

yy  L xx  L K yy  . 2 y yy J 


M = D (1-v)  h 
xy  xy 

where 

E = Et/(l-v2);  D = Et3/12(l-v2);  X =A  (l-v2)/tX 

XX  XX  X X 

p 

X = A (l-v2)/ti.  ; p = El  /VI  ; and  p = El  /DZ  . 
yy  v'  y Kxx  xc  x Kyy  yc  y 

From  Eqs . (2)  one  may  derive  the  following  expressions  for  the 

reference  surface  strains 

e = a.N  + a„N  + a_x  + a.x. 
xx  1 xx  2 yy  3 xx  4 yy 


e = a„N  + b N + b n + b h 
yy  2 xx  2 yy  3 xx  4 yy 


(3) 


where 


e - hy  ~ N /( l-v)E 
xy  xy  xy  xx 


al  = (1+Xyy)/Q'Exx  ; a2  = _v/aExx  ; *3  = ( 1+Xyy  )exXxx/c> 

P P 


a.  = -ve  A /a\  b.  = (1+X  )«E  ; b = -ve  X la  (4) 

4 y yy  2 xx  xx  3 x xx 
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b - ( 1+A  )e  \ /a',  a ~ [ 1+X  )(1+A  ) -v^] 

4 xx  y yy  xx  yy 


By  employing  the  principle  of  the  stationary  value  of  the  total  potential 

one  can  derive  the  following  equilibrium  equations 

N + N =0 
xx,x  xy,y 


N + N =0 
xy,x  yy ,y 


N 


M + 2M  + M = + [N  (w,  +w°  )],  + [N  (w,  +w"  )], 

xx, xx  xy.xy  yy.yy  R yy  y y y xy  x x y 

+ [n  (w,  +w°  )],  + [N  (w,  +w°  )] , 

xx'  ’x  xx  xy  y y x 


By  introducing  the  Airy  stress  function,  as  N = -N  + F,  , N = F, 

° J xx  xx  yy  yy  xx 

anti  N = -F,  where  N is  the  level  of  the  applied  uniform  axial  com- 
xy  xy  xx 

pression,  the  first  two  of  Eqs.  (5)  are  identically  satisfied. 

Next,  by  eliminating  u and  v from  the  first  three  of  Eqs.  (1),  employ 
ing  Eqs.  (3),  the  Airy  stress  function  and  the  last  three  of  Eqs.  (1)  one 
can  derive  the  compatibility  equation  in  terms  of  the  Airy  stress  function 
F and  the  radial  displacement,  w.  If  one  expresses  the  third  of  Eq . (5) 
in  terms  of  F and  w,  the  governing  equations  consist  of  two  coupled  par- 
tial differential  equations  in  F and  w.  These  arc: 

Equ i 1 ibr ium 

DL.  [w]  - L [F]  - F,  /R  + N (w,  + w°  ) -L[F,w  + w°]  = 0 (6) 

n q xx  xx  xx  xx 


Compat ibi 1 ity 


Ld[F]  + L [w]  + '.>L[w,w  + 2w  ] + w,xx/R  = 0 


(7) 


where  , and  L(  are  differential  operators  defined  by  L , 
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Lg[sl  8ns’xxxx  + 2g12s,xxyy  + S22S ’yyyy 


d..  = (1+X  )/aE 
11  XX  XX 


di2  ■ [<1+\«)<1+V  • v:l/‘'(1-v)Exx 


- ( 1+2.  )/oE 
22  yy  xx 


10  e2X  (1+X  -v2) 

12  x xx  yy 


h.-l  + p + ^ ^ ^ 

11  xx  j.2  a 

ve  e X X 

h = i + 11  * y .**  yy 

"12  2 Q 

t 2 2 

eX  (1  + X - v ) 

h = 1 + d + - ^ Si 

h22  1 pyy  t2  a 


a. . = - ve  X /a 

411  X XX 


q = (1  + X )eX  + (1+X  )e  X ]/(2o) 

4 12  yy  x xx  xx  y yy 


q22  = ‘ vey  V“ 

and  L is  a differential  operator  defined  by 

LrS,T]=S,  T,  - 2S,  T,  + S , T,  (8b) 

L ’ J ’xx  yy  ’xy  xy  yy  xx 

The  total  potential  expression,  in  terms  of  the  Airy  stress  function  and 


the  radial  displacement,  is  given  below 


UT  ■ 2T-  I.  <£ir’yy  + 82F'xx  + 83F,xxF,yy  + V'xy)dA 


'xx  A 
P 


(10) 


+ \ JA  (Ofjw^yy  + a2w^xx  + c3v,xxw)yy  + o'4w^xy)dA 

N P,  2 

- r (23.F,  + ?,F,  )dA  + ttRL  N - N 2nRLc.„ 

2E  J,  1 yy  3 xx  E xx  xx  AV 

xx  A xx 


where  e (average  end  shortening)  is  given  by 

A V 


e. . - - u,  dA/2nRL 
AV  J . x 
A 


and 


21  d22Exx  ’ 32  dllExx  ’ 5 3 " 2v/cv;  24  2/(1  ‘ v) 

P P 

]9  e e X X 

Oi  = h22;  a2  = hu;  o3  = 2^  1 + ±|  — ^x-x-^-v ; a4  = 2(1  - v)  (11) 

Similarly,  the  expressions  for  the  average  end  shortening  and  "unit 
end  shortening"  at  y = 0 are  given  by 


1 


2ttR  L 


e,„  = a.N  - 7 — La.F,  + a„F,  + a.w,  + a.w, 

AV  l xx  2ttRL  • •>  1 yy  2 xx  3 xx  4 


yy 


w,  (w,  + 2w',  )1dxdy 

2 x x x 


(12a) 


e = a.N  - - , a F,  + a F,  + aw,  + a.w, 

1 xx  L • L 1 yy  2 xx  3 xx  4 yy 


- 2 w>x  (w’x  + 2w’x}  dx 

y=o 


(12b) 


Note  that  e measures  the  amount  of  end  shortening  per  unit  of  cylinder 
length,  L.  The  associated  boundary  conditions  are  either  kinematic  or 
natural  (but  not  both)  except  for  the  direction  associated  with  the  length 
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of  the  cylinder  (x)  in  which  case  the  displacement  component  u is  free 


and  the  stress  resultant  N must  equal  to  the  applied  stress  resultant, 

xx 

-N  . Thus  at  a boundary  characterized  by  x = 0 or  x = L the  boundary 
xx 

conditions  are: 

Either  or 


in-plane 


N = F , - N =-N 

xx  yy  xx  xx 


u = constant 


N =0 
xy 


v = 0 


transverse 


M =0 
xx 


w . =0 

X 


* 

Q = 0 

x 


w = 0 


The  expressions  for  the  moment  resultant,  M and  the  effective  transverse 
shear,  , are 


,v,w,  + Y-%w>  + v0(F,  - N ) + Y/Fj 

Y1  xx  2 ?yy  Y3  yy  xx  r4  xx 


o"  = (F,  - N )(w,  + w”  ) + F,  (w,  + w°  ) - M - 2M 

<x  u ’yy  xxA  ’x  ’x'  ’xyv  ’y  ’y  xx,x  xy,y 


where 


Yl  = Dhn;  y2  = D — ; \3  ~ - ay>  Y4  = -b: 


The  general  computer  program  is  written  for  the  following  end  conditions 
(SSi,  CC  i ,FFi  , i = 1,2, 3, 4) 


I 


SS  : 

: w = 

M =0 

XX 

1. 

F,  = F,  =0 

xy  yy 

cc 

: w = 

w,  =0 

X 

2. 

F,  =0;  u = C 

xy 

FF: 

* 

: Qx 

= M =0 

XX 

3. 

< 

II 

II 

o 

(15) 

4. 

< 

II 

o 

c 

II 

o 

where  C is  a constant. 

The 

conditions  in  u and 

v can  be  expressed  in  terms  of  w and  F 

as 

in 

Ref. 

14.  For  example,  the 

condition  u = C 

in  SS  -2  can  be  replaced  by 

a condition  expressed  solely  in  terms  of  w,  w , F and  their  gradients. 

This  is  accomplished  by  the  following  procedure: 

This  boundary  condition,  SS=2,  at  x = 0 or  L given  w = 0;  F,  = 0, 

M =0  and  u = C.  The  first  two  are  in  terms  of  w and  F.  The  third  one, 
xx 

= 0,  from  the  first  of  Eqs . (14)  is  expressed  in  terms  of  w,F,  and 
their  gradients,  ter  the  last  one,  one  notes  that  [see  Eqs.  (1)  and  (3)] 

1 


e = ~ Lu.  + v , 
xy  2 y x 


+ w,  w,  + w,  w°,  + w,  w°  ] = -F,  / 

x y y ’x  x y ’xy  (1-v 


)E 


P 

(16) 


since  F,^  = 0,  w,^  = 0 because  w(0,y)  = 0,  and  u,  = 0 because  u(0,y)  = C, 
Eq . (16)  becomes 


v,  + w,  w, 
x x y 


(17) 


Similarly,  from  Eqs.  (1)  and  (3)  one  may  write 


e = v,  + s[w,  (w,  + 2w°  )]  - 'i 

yy  ’y  2L  ’yv  ’y  V J R 


= a„N  + b„N  + b k + b.  k 
2 xx  2 yy  3 xx  4 yy 


(18) 


This  equation,  Eq.  (18),  is  valid  at  any  point  along  the  shell,  thero- 
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fore  differentiation  with  respect  to  x does  not  violate  its  validity. 
If  this  is  done  and  if  the  N's  and  k's  are  expressed  in  terms  of  w,  F, 


gradients 

+ l [w, 

2 xy 

, one  may  write 
(w , + 2w°  ) + w , 

(w. 

+ 2w' 

y y y 

’xy 

= aoFj 

+ b„F , + b w, 

+ b. 

w. 

2 yyx 

2 ’xxx  3 xxx 

4 

yyx 

(19) 


Evaluation  of  Eq . (19)  at  x = 0 or  L,  and  use  of  the  fact  that 
w>,,  (0,y)  = 0 yields 

v,  + w,  w°,  - w,  / = a„F,  4-  b F,  + b w,  + b w,  (20) 

yx  xy  y x R 2 ’yyx  2 ’xxy  3 xxx  4 yyx 


Differentiation  of  Eq.  (17)  with  respect  to  y,  yields 


. o . o 

v,  + w , w,  + w,  w , =0 

xy  xy  y x yy 


(21) 


Substitution  of  F.q.  (21)  into  Eq.  (20)  yields  a boundary  condition  equi- 
valent to  u = C,  or 


b F,  + b w,  + b.w,  + w,  (^  + w"  ) = 0 
2 ’xxx  3 ’xxx  4 yyx  x R ’ 


yy 


(22) 


Similar  steps  may  be  followed  to  express  all  possible  boundary  conditions 
in  terms  of  w,  F,  and  their  gradients.  In  order  to  save  space  only  the 
final  expression  for  all  possible  boundary  conditions,  Eqs . (15),  are 
given  below,  which  have  been  incorporated  into  the  computer  program  (see 
Appendix  B) 


SS-1  w = v.w,  + v.F,  — F,  - F,  - 0 
'1  xx  ’4  xx  xy  yy 


' F 

3\  : 

, - N ) 

yy  xx 

+ y4f, 

+ h.w, 

4-  w, 

xxx 

4 yyx 

X 

xy 
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SS-3 


w = y,w,  + y.F,  = F,  = b F,  + b w,  = 0 
'1  xx  '4  xx  yy  2 ’xx  3 xx 


SS-4  w = v,w,  + v^(F  - N ) + v,F,  - a„(F,  - N ) + b„F, 

T1  xx  T3  yy  xx  T4  xx  2 yy  xx  2 


+ b w , - 0 

3 xx 


a + 2/  F,  + b F,  + b w,  + t w, 

2 (l-v)E  J xyy  2 xxx  3 xxx  4 xyy 

XXP 

' + w,  (]  + w°  ) = 0 
x R yy 


(23a) 


CC-1  w=w,  = F , = F , *0 

x xy  yy 


CC-2  w = w,  = F,  = 0 
x xy 


b F,  + b w,  = 0 
2 xxx  3 xxx 


CC-3  w = w,  = F,  =0 
x ’yy 


b.F,  + b w,  = 0 
2 xx  3 xx 


CC-4  w = w,  = a „!  F, — N ■ + b„F,  + b„w,  = 0 
x 2'  yy  xxy  2 xx  3 xx 

a0  + 2 / ( l-v)E  F , + b_F , + b.w,  = 0 

2 xx  i yyx  2 xxx  3 xxx 

P 


(23b) 


S imilar ly  the  FF-1  condition  and  the  symmetry  and  antisymmetry  conditions 
at  x = L/2  are 


FF-1  v,w,  + v^w,  + y,F,  = F,  - F,  - 0 
' 1 xx  T2  yy  T4  ’xx  yy  xy 


/.F,  + Yiw»  + Yt  + 2D(1-v)  w,  + N (w,  - w,  ) = 0 

4 xxx  ’1  xxx  ^'2  _i  xyy  xx  x ’x 

(23c) 


Symmetry  ^w,x  = ^ = ^xy  = u = 0) 
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= -Y.w, 

' 1 xxx 

+ 

Y/F,  - (F, 

T4  xxx  yy 

= b F, 

+ 

b w , 

o 

+ w,  w. 

2 ’xxx 

3 xxx 

’x  *' 

Antisymmetry  (w  = M = v = F.  = 01 
‘ ^ xx  yy 


w - v w , + y . F , =0 

'1  xx  '4  xx 


(23d) 


F,  - b F,  + b w,  = 0 
yy  2 ’xx  3 ’xx 


The  problem,  as  formulated  herein,  is  to  find  the  limit  point  which 
represents  the  buckling  load  for  the  imperfect  configuration.  This  implies 
to  solve  the  field  equations,  Eqs . (6)  and  (7),  subject  to  the  proper 
boundary  conditions  for  a given  imperfection  and  level  of  the  applied  load 
(small  initially),  -N  , and  thus  obtain  the  corresponding  amount  of  "unit 
end  shortening",  Eq . (12b).  By  plotting  N versus  e one  can  obtain  the 
limit  point  (theoretically). 
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CHAPTER  III 


SOLUTION  METHODOLOGY;  NONLINEAR  BUCKLING  ANALYSIS 
By  employing  the  von  Karman-Donnell  kinematic  relations  the  field 
equations  consist  of  two  coupled,  nonlinear,  partial  differential  equa- 
tions in  terms  of  the  transverse  displacement,  w,  and  the  Airy  stress 
function,  F.  The  procedure  employed  herein  for  accomplishing  a solution 
is  basically  similar  to  that  of  Refs.  13  and  14.  The  system  of  partial 
differential  equations  is  reduced  to  a system  of  ordinary  differential 
equations  by  using  a separated  solution  (Fourier  series)  of  the  following 
form  (see  Refs.  13  and  14). 

K 

w(x,y)  = £ W.  (x)  cos  pp 
i=0  1 R 

2k 

F(x,y)  =.§0  f i (x)  cos  ijP  (24) 

K 

w°(x,y)  = Z W°  (x)  cos 

i=0  1 R 

Note  that  W°(x)  denotes  the  known  coefficient  of  the  ith  component  of  the 
geometric  imperfection,  and  n is  the  parameter  associated  with  the  number 
of  full  waves  around  the  circumference. 

By  substituting  Eqs . (24)  into  Eq . (7),  employing  trigonometric 
identities  of  double  Fourier  series  as  in  Ref.  19  involving  products  and 
the  orthogonality  of  the  trigonometric  functions,  the  compatibility  equa- 
tion becomes 
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for  i = 0 


£o  - 17-  - - »°/R  + 4 i J2(«  + Wj)», 

11  4R  j=1  J J J 


(25a) 


for  i = 1,2, . ■ ,2K 


dllf-  - 2' 


i - At)  di2fi  + vt'J  d22fi 


in/4 


""  'in\2  ..  " . /in4 


+ 6iLqnwi  " 2It;  c‘i2wi  + q22wi + Wi/R 


n2  K 0 

n , , /f  . .,2 


772  £ h+J<“l+J  + 2“°+J> 


(25b) 


+ (2-4i)(t-J,2V-)1<“U-)l  + 2Wfi-j|)J"j 

n on  9 M 

+ f 6 . . (W  + 2WT,  .)  + (2  - T)  . . ) 6 1 . I (W|  . 1 o"  ,2 

i+J  i+J  i+J  'j-i'  | i-j|V  I i"  J I + 2W|  .|  )]jZW. 

+ 2Ul+j)t.  + ,(«:  + . + 2U°^.)  - 'l1.J|l-j|6|1.J|(W|'...|  +w|’j.J|)]jwjj  -0 


where 


6f  = 


V1 


l > K 


^ K 


- 1 a < 0 

Tlx  = ^ 0 X = 0 

1 X > 0 


and 


( ) = — 
V ’ dx 


Next,  if  Eqs . (24)  are  substituted  into  the  equilibrium  equation,  Eq.  (6), 
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iny 

and  the  Galerkin  procedure  is  employed  (cos  -tt-  is  the  weighting  function 

K 

i = 1,2,...K),  the  vanishing  of  the  (K+l)  Galerkin  integrals  leads  to  the 
following  system  of  (K+l)  ordinary  differential  equations 
for  i = 0 


im  ■)  n 7-0 

Wo  [Dhu  + qn  /dn]  + Wo[2qil/R-dn]  + Wo[l/R  * d^]  + 


2 K 


■2—  Z j2i-ii[(w.  + 2W°)W'.'  + (WV  + 2W°")W. 
4R2  j=l  ldll  J J J J 


+ 2(W'.  + 2W0,)W'.]  + r (W.  + 2W°)W.]  - 2[(W.  + W°)fV 
J J J Rdn  J J J J J J 


+ CM"  + W°")f . + 2(W!  + W°')f  !])  = 0 
J J J J J J J 


(26a) 


for  i = i ,2  , . . . ,K 


»>n»r  - z( T 2 h12ui'  + (r)4  h22“lj 


L’ufr  - 2i  f)  ’i2f;' + (tv  ’22fi + f,i,/R 


( in\^ 


+ N (WV  + w°")  - W'  (ij?)  ! (w  + w°)  - w f-±-(  ^ (W 

XX  i l 0L-dll  ^ K J i 1 oLRd^v  R / v l 


4 .2  K 

+ W°)  + -2-  (W.  + W°)  E j (W.  + 2W  )W. 

l , „4  d , l i J J J 


4R  U11 


j=l 


2 2K 


f o 0 

+ 2T2  j5li[(i+j)  6i+3^+.i+Wi+j) 
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+ (2  - T]2  .)  ( i-  j ) 2 6 . . (W,  . + W°,  . 

J-i  |i-j|  I 1-J|  I 1-JI  J 


+ [6.^.(W''  . + W°")  + (2  - T)2  -)6|  . .|  +wfV 

L l+jV  l+J  1+j'  'j-i  | x- j | |i-j|  |l-j|  j 


+ 2[  (i+j)6i+J  (w:+j  +w°;.) 


T'i-jli"jl6|i-j|(Wfi-j|  +W|i-j|)]  jfj)  = 0 


(26b) 


For  a given  value  of  the  applied  load,  N , and  imperfection,  Eqs . (25) 
and  (26)  represent  a system  of  (3K  + 2)  coupled  nonlinear  differential 

equations  in  (3K  + 2)  unknowns,  f^  i = 0,1,2,... .2K  and  i = 0,1,2 ,K. 

These  equations  denote  equilibrium  and  compatibility  conditions.  Similarly, 
the  expressions  for  the  total  potential,  UT,  average  end  shortening,  eAV> 
and  "unit  end  shortening",  e,  become 
L 


U = nR 
T Jo 


/ , .po  _ W 2 ~ „ ,2 

( -1—  ! - —■  - q-.W"  + 1 (W  . + 2W?)W  . 

\ E 1.2  L R M11  o ._2  l i iJ 

\ xxp  dll  4R  5=1 


2K 


+ 2 I lBl(f)4  fi  + hCf  - »&f  fifi  + »*(£)’  fl2]} 

i = l 


+ + i 


K 

1=1 


2 


2 


ai(ir)  wi+wf  - «3(t)  w±wi + ^(t)  w±2]) 
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N P, 
xx  3 

d.  ,E 
1 1 xxp 


W 2 K 

f'  qHWo+^  I i2(W.  +2W 
4R  i=l 


°>"d) 


dx 


P.ttRL  _2 

+ N - N 2rrR  LcA17 

E xx  xx  AV 

xxp 


(27) 


K 

V .2. 


, L a -W  2 

e = a.N  + f [ -i  -r=-  | — + q . .W  " - — ) i"  (W  . + 2WV)W  . 

AV  1 xx  L I d17  L R nll  o .2  Z_i  l l l 

oil  4R 


i=l 


K 


- a„W " + W 1 (W  ' + W°')  + y w : (W  : + wV')|-  dx 
3o  2 o o o 4 _.  l i l J 

i=l 


(28a) 


, L / a,  -W 
1 r>  / 2 o 


K 


e = a.N  +71  ( - — — + q11W" i (W  . + 2W  . )W  . 

lxx  L Jo\  dn  LR  ^11  o ,„2  1 ’ 1 


J11 


" ^ -V  ■ 

1 1 1. 


4R  I=i 


2K 


r /in' 


• \ 2 


K 


K 


+ ai[ir)  f-  ■ a9f"  - wv  + a w. 

IV  R / l 2 iJ  3 /_  i 4 V R / i 
i=l  i=0  i=l 


K K 

W ‘ 

i=0  i=0 


: (w!  + 2W°’)  \ dx 


(28b) 


Finally,  the  appropriate  boundary  conditions  are  expressed  in  terms 

of  W , W°,  and  f..  Note  that,  because  of  the  character  of  the  nonlinear 
l i l 

differential  field  equations,  Eqs.  (6)  and  (7),  by  setting  n = 0 one 
obtains  the  linearized  version  of  equilibrium  and  compatibility.  Further- 
more, it  is  easily  seen  from  Eqs.  (24)  that  n = 1 includes  the  axisymmetric 
mode  since  the  summation  on  _i  starts  from  zero.  In  addition,  it  is  seen 
from  the  Fourier  series  representation  of  the  imperfection  that  this 
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expression  is  suitable  for  the  case  when  the  imperfection  is  of  the  same 
shape  as  the  buckling  mode, as  well  as  for  any  arbitrary  symmetric  (with 
respect  to  y)  imperfection  shape.  In  this  latter  case,  in  order  to  obtain 
a solution  it  is  necessary  to  let  n = 1 (in  the  series  representation 
for  W°)  and  take  K large  enough  for  an  accurate  representation  of  the 
imperfection  and  for  achieving  a convergent  solution.  By  employing  Eq . (24) 
the  expressions  for  the  various  boundary,  symmetry  and  antisymmetry  condi- 
tions become,  Eqs . (23), 


SS-1  W = W"  = 0 

o o 


W i - y^W V + y,  f"  - 0 i - 1,2,  ...,K 


f!  = f.  = 0 

l l 


i = 1,2,. ..,2K 


W = W"  - y^N  = 0 
o o ' 3 xx 


“i  ’ - Y3l*«  + (f)  ft]  + Vi  ’ 0 1 ' l"2 K 


f.'  = 0 
1 


i = 1,2,  ...,2K 


vi" + Vi"'  - (f)  v;+iK 


> 1 (i+j)V\  + (1  - "|*  . + 71  ) ( i ~ j ) 2W ? . j iw!  = 0 (29a) 

2r2  1=0  L i+J  J 1 ' 1 J'  J 

i = 1 ,2 , . . . ,2K 


w = w"  = o 

o o 


w.  = Y,W''  + Y,fV  = 0 
t Y1  l 4 l 


i = 1 ,2 , . . . ,K 


I 


f.  = b„f7  + b WV  = o 

l 2 l 3 i 


i = 1,2,... ,2K 


SS-4 


W = W"  = 0 
o o 


= v W"  - vJ  I ^r)  f.  + N | + y.t'!  _0  i 1»2,...,K 

i ’1  i yy_\  R ) 1 xxJ  T4  l 


Tn  + £.  ; + b0f'.'  + b.WV  = 0 ; 

a2L  xx  \ R / iJ  2 l 3 i 

2 


in\ 


^ + 2/(l-v)ExxpJ  v r 


f:  + b2f-  + b3w-  - b4^Tj  W! 


= 1,2,. ..,2K 


/ in\  ^ 

^ w:  + 


si  ^ [<i+j)V+j  + (1  - -n^-i 

2R  :_n 


j=0 


2 + T1 . )W? . . I ! w ! = 0 

'l  1-J  J J 


i = 1,2,...  ,2K 


CC-1 


w . = w ! =0 
1 1 


f!  = f.  = 0 
1 1 


i — 0,1 ,2  , . . . ,K 
i = 1,2,  ...,2K 


CC-2 


w . = w ! =0 
1 1 


i ~ 0,1 ,2  , . . . ,K 


CC-3 


n “Vi  +b3Wi  = 0 


W.  = w!  =0 
1 1 


f.  = b.f'.'  + b W'.'  = 0 
1 2 1 3 1 


, i = 1,2 2K 

i = 0,1 K 

, i = 1 ,2 , . . . ,2K 


(29b) 


CC-4 


W = W : =0 

1 1 


i “ 0,1 ,2  , . . . ,K 


22 


I * 


a , N + {—)  f . + b f'.'  + b W"  = 0 

21.  xx  \ R / 1 j 2 l 3 l 


i = 1,2, ...,2K 


'*2  + 2/<1-V)Exxpj  (ffh  * W + b-WV  - 0 , 

3 l 


i = 1,2,...,2K 


V r .2, 


- VHi/diU  - V4/Rdii  * -J—  - s f“j  + “j’  “j  ■ 0 


4R  du  j = l 


V,  - V/ q i i /d , i ! + w'  , N - y. /Rd . . + N W° ' 

o l.1  T 4 n 1 1 11  _ o !_  xx  T4  11J  xx  o 


V4n 


K 


4R2d 


11  j = l 


j (W J + 2W°)JW^  + (Wj  + 2W°')vr  j = 0 


^ih;  - (f)  Y2wi + v 'i  - ° 


i 1 ,2  , . . . ,K 


(29c) 


Vi"  +^iwi' 


N 2 

Y,  + 2D  ( 1-v  ) ; W'.  + N (W ! + W°')  = 0 

.2  J \ R / j xx  l l 


i = 1 ,2 , . . . ,K 


f!  = f.  = 0 

l l 


i = 1,2,... ,2K 


Symmetry 


W'  = 0 
o 


C^l  ' Vll^lP  + “2 


4R  d 


11  j=l 


j2'  (W.  + 2W°)W! 
L i j j 
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2 

+ (Wl  + 2W°')W  ' + 2—  ; I 6.  ,W°'  + (1-T12  . + 71 . )W? ! 6,  . j2f  =0 

J J J 2R2  L > 1+J  1+J  1 u-jI  U-JI-  J 


2K 


W!  = Yl«r  + y4f-  + ^ I + (1  - ^ + V 6|t.jlwj,;.jl]jf1 

j=l 


j=l 

2K 


i = 1 ,2 , . . . ,K 


(29d) 


f!  = b f V 1 + b.W'." 
i 2 l 3 i 


2R 


2 L L 

j=i 


W°,  . + (1  - Tl2  . + T1 . )W? ! . = 0 

i+J  J-i  'i  |i-j|J  J 


i = 1,2,.. .,2K 


Antisymmetry : same  as  SS-3  (29e) 

The  solution  procedure  employed  is  described  below. 

A generalization  of  Newton's  method  (Refs.  20  and  21),  applicable  to 
differential  equations,  is  employed  to  reduce  the  nonlinear  field  equations, 
Eqs . (25)  and  (26),  and  appropriate  boundary  conditions  to  a sequence  of 
linear  systems.  In  this  method,  the  iteration  equations  are  derived  by 
assuming  that  the  solution  is  achieved  by  a small  correction  to  an 
approximate  solution  (initially  taken  as  the  linear  solution).  These 
small  corrections  are  obtained  through  the  solution  of  the  linearized 
(with  respect  to  the  corrections)  differential  equations. 

The  linearized  differential  equations  are  written  in  matrix  form  as 
follows : 

Field  equations 

[R]  fZ"}  + [S]  {Z ’ } + [T]  fz}  = fg}  (30) 
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Boundary  Conditions 

[sl  fz'3  + [T]  fz]  = fg] 


where  fzj  is  the  vector  of  the  6K  + 2 unknowns. 


(31) 


fz}1  = fWo,W1>...WK,f1,f2,...f2K,  f'^f^.-.f^]  (32) 

Note  that  f has  been  eliminated  in  a manner  similar  to  that  of 


o 

Refs.  13  and  14. 

These  ordinary  differential  equations  are  cast  into  the  form  of  finite 
difference  equations,  and  the  system  of  ordinary  equations,  Eqs . (30)  and 
(31),  are  changed  into  a system  of  linear  algebraic  equations.  The  usual 
central  difference  formula  is  used  at  all  mesh  points,  i.e. 


= (Z 


X+l 


- Z 


Z-l 


) /2A 


(33) 


Z'[  (Z£- 1 ‘ 2Zl  + Z£+l)/A 

Note  that  the  second  derivatives  in  W.  and  f.  are  taken  as  independent 

it 

elements  of  the  vector  of  the  unknowns,  therefore  the  second  of  Eqs.  (33) 
applied  only  to  fourth  derivatives  of  VL  qnd  f^.  By  using  one  fictitious 
point  on  each  side  of  the  cylinder  ends  one  obtains  a system  of  (6K  + 2)  X 
(NP  + 2)  difference  equations  (NP  - number  of  mesh  points). 

These  equations  are: 


Cci:  fzj  + [B1]  fzj  + [a l1  fz2]  = ^ 

[cxD  fz£_L)  + [b£]  fz,}  + [a,]  fz£+1]  = g,;  x = 1,2,...,NP 

1 Sp"  f ZNP- 1 ^ + l\p^  fZNP^  + ^NP^  'ZNP+1^  ” ^IP 
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whore  fZQ]  and  [Z^p^j  are  unknown  vectors  at  the  iwo  fictitious  points, 
and  the  matrices  in  Eqs . (34)  are  given  by 

L^r  = ' 2A  ^l-1  J ^NpJ  = ‘ 2A  ^NP-* 


[Bx]  = [?{!  ; [§Np]  = [tnp3 

'-^1J  " 2A  ’ ^NP-^  = 2A  ^NP' 


This  system  is  solved  by  the  special  algorithm  which  is  reported  in  Ref.  22. 

When  the  load  parameter  is  at  a limit  point  a unique  solution  does 
not  exist  (the  system  becomes  singular)  and  thus  the  solution  does  not 
converge.  Therefore,  the  solution  procedure  goes  as  follows:  first, 
the  system  of  equations  is  solved  for  a small  level  of  the  applied  load 
(say  20%  of  the  classical  buckling  load),  then  a multiple  of  this  solution 
is  used  for  a small  increase  in  the  load  parameter  until  the  process  fails 
to  converge.  The  load  level  at  which  the  solution  fails  to  converge  is 
taken  to  be  the  critical  load.  Note  that,  when  approaching  the  limit 
load,  if  the  increment  in  the  load  value  is  large  enough  so  as  to  place 
the  systems  at  equilibrium  far  beyond  the  limit  point,  the  system  in  some 
cases,  does  converge.  In  this  case,  since  the  interest  is  in  the  limit 
point  value  only,  one  can  check  the  sign  of  the  determinant  of  the  coeffi- 
cients of  the  unknown  vector.  If  the  sign  changes,  by  taking  a large 
increment  in  the  load  level,  one  must  decrease  the  increment  and  proceed 
with  the  solution.  Large  increments  are  used  in  the  procedure  in  order 
to  save  computer  time.  Because  of  the  use  of  large  increments  and  since, 
in  some  cases,  the  solution  converged  at  both  consecutive  steps  the 
criterion  of  the  change  in  the  determinant  sign  is  employed  to  establish 
the  existence  of  a critical  point  within  the  range  of  these  consecutive 
steps.  At  each  level  of  the  load  for  which  the  sytem  is  solved,  the 
value  of  the  number  of  full  circumferential  waves  is  needed.  Different 
values  of  n are  used  to  obtain  a solution,  and  the  one  that  minimizes 
the  total  potential,  Eq.  (27),  is  taken  as  the  correct  one  (see  Refs.  13 
and  14).  Numerical  integration  is  used  to  find  the  total  potential. 
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The  number  of  n values  to  be  tried  at  every  increment  of  the  load  is 

small,  since  the  circumferential  mode  does  not  vary  significantly  with 

small  increases  in  the  load,  N . The  end  shortenings  at  each  level 

xx 

of  the  applied  load  are  also  computed  through  numerical  integration. 


Chapter  IV 


APPLICATIONS  AND  DISCUSSION 

The  mathematical  formulation  and  the  method  of  solution  for  the 
buckling  analysis  of  imperfect,  thin,  circular,  stiffened  cylindrical 
shells  under  uniform  axial  compression  is  presented  in  Chapters  II  and  III. 
The  methodology  is  demonstrated  through  a number  of  illustrative  examples. 
Numerical  solutions  are  obtained  by  employing  the  Georgia  Tech  high  speed 
digital  computer  CDC-CYBER  70,  Model  74-28.  A general  program  is  written 
(see  Appendix  B)  which  includes  the  following  desirable  features:  (a) 

it  is  applicable  to  stiffened  (in  either  or  both  directions)  geometries 
as  well  as  unstiffened;  (b)  it  accomodates  all  possible  boundary  condi- 
tions (SS,CC,FF,  etc.),  and  it  can  easily  be  modified  to  accomodate  elastic 
end  restraints;  (c)  the  number  of  Fourier  terms  (K)  can  be  as  large  as 
desired.  The  same  holds  true  for  the  number  of  points  (NP)  in  the 
finite  difference  scheme;  (d ) the  geometric  imperfection  can  be  axisym- 
metric  as  well  as  an  arbitrary  symmetric  (w.r.t.  y)  one.  The  program 
can  easily  be  modified  to  include  other  destabilizing  loading  conditions 
such  as  pressure  and  torsion. 

Although  the  program  is  highly  dimensional  because  of  the  number  of 
Fourier  terms,  number  of  points  and  number  of  required  iterations,  the 
solution  is  obtained  with  reasonable  CPU  time.  For  example,  by  using 
K = 1 (one-term)  and  65  points  (536  unknowns)  it  requires  four  seconds  to 
complete  one  iteration;  for  the  same  but  K = 2 it  requires  15  seconds. 

For  a convergent  solution  to  be  obtained  at  low  levels  of  the  applied 
load  two  iterations  are  sufficient.  At  load  levels  approaching  the  limit 
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point  six  iterations  are  needed  (convergence:  percent  difference  <10  ). 

The  numerical  results  for  all  the  illustrative  examples  are  presented 

in  tabular  form  in  Table  1.  A number  of  these  examples  are  taken  from 

the  open  literature  in  order  to  check  the  present  solution.  In  addition, 

new  results  are  generated  and  the  discussion  of  both  is  given  below.  In 

this  table,  for  each  example  considered,  the  buckling  load  (see  columns 

of  N and  N /N  ) is  bracketed  between  two  numbers  (denoting  the 

xx  xx  XX  „ ve. 

cr  cr  c l 

desired  accuracy).  The  first  number  denotes  the  highest  level  of  the  load 
for  which  a convergent  solution  is  obtained  and  the  second  number  a level 
higher  than  the  limit  point  (according  to  criterion  discussed  herein). 

In  all  examples,  the  imperfection  is  taken  to  be  symmetric  with  respect 
to  x = L/2  and  therefore  the  response  is  taken  to  be  symmetric.  Thus 
only  half  of  the  cylinder  is  analyzed  by  employing  the  appropriate  symme- 
tric conditions  at  x = L/2.  The  examples  can  broadly  be  classified  in 
one  of  the  following  four  categories:  (a)  unstiffened  (Examples  1-5); 

(b)  stringer-stiffened  (Examples  6-9);  (c)  ring-stiffened  (Examples  10 

and  11);  and  (d)  ring-  and  stringer-stiffened  (Examples  12-21).  In  each 
of  the  four  categories,  at  least  one  case  was  calculated  for  two  different 
truncated  Fourier  series  (K  = 1 and  K = 2)  and  for  two  different  number 
of  points  (NP  = 35  and  NP  = 65)  in  order  to  check  the  effect  of  these 
two  parameters  on  the  convergence  of  the  solution. 

(a)  Unstiffened  (Examples  1-5).  This  geometry  is  taken  from  Ref.  16. 

Only  Example  3 is  reported  in  Ref.  16  and  the  results  are  in  very  good 
agreement.  Examples  1,  2,  4 and  5 are  considered  in  order  to  assess  the 


effect  of  boundary  conditions  for  this  type  of  an  imperfection  (see  Table 
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1).  It  is  seen  from  the  results  that  the  effect  of  in-plane  boundary 
conditions  is  significant  for  the  simply  supported  case  and  insignificant 
for  the  clamped  case.  In  Example  2,  because  of  the  u = C in-plane  boun- 
dary condition,  the  criterion  used  for  finding  N is  that  the  deter- 

cr 

minant  (see  Chapter  3)  goes  to  zero.  It  can  be  seen  from  Fig.  2 that 
as  the  sign  of  the  determinant  changes  the  radial  mode  of  deformation 
changes.  The  calculations  for  this  case  were  based  on  K = 1 and  65  grid 
points  in  the  axial  direction  (for  half  the  cylinder  length). 

(b)  Stringer-stiffened  (Examples  6~9).  The  geometry  for  these  examples 
is  taken  from  Ref.  9 and  referred  to,  in  it,  as  heavy  stringers.  The 
present  results  are  in  very  good  agreement  with  those  of  Ref.  9.  Examples 
6 and  7 correspond  to  Z = 95.4  with  external  stringers  and  the  imperfec- 
tion shape  is  virtually  axisymmetric  in  Example  6,  and  symmetric  in 
Example  7.  The  axisymmetric  imperfection  yields  greater  reduction  than 
the  symmetric  one.  Note  that,  for  this  case,  the  perfect  geometry 
buckles  axisymmetr ically . Examples  8 and  9 correspond  to  Z = 394  with 
external  and  internal  stringers  respectively  and  an  imperfection  which 
is  primarily  axisymmetric.  From  these  examples  one  can  conclude  (as  in 
Ref.  9)  that  imperfection  sensitivity  is  greatly  affected  by  the  curva- 
ture parameter  Z,  and  that  externally  stiffened  configurations  are  much 
more  sensitive  to  geometric  imperfections  than  internally  stiffened  ones. 
These  cases  were  analyzed  with  35  and  65  grid  points  in  the  axial  direc- 
tion for  half  the  cylinder  length.  From  these  computations  one  can  see 
(for  a typical  case,  see  Tables  2 and  3),  that  the  critical  load  and  the 
response  (w,F)  for  each  n are  almost  the  same  regardless  of  the  number 


31 


r 


1 

of  points. 

(c)  Ring-stiffened  (Examples  10  and  11).  This  geometry  is  also  taken 
from  Ref.  9,  and  it  corresponds  to  light  ring-stiffened  cylinders  with  Z 
= 394  and  external  and  internal  positioning  of  the  rings  respectively. 

Although  the  numerical  results  of  the  present  analysis  are  in  good  agree- 
ment with  those  of  Ref.  9,  the  conclusion  concerning  the  effect  of  ring 
positioning  on  the  sensitivity  is  reversed.  According  to  the  present 
results  internal  rings  make  the  overall  configuration  more  sensitive  than 
external  rings. 

(d)  Ring  and  stringer  stiffened  (Examples  12-21).  These  examples  are 

chosen  to  demonstrate  the  methodology  for  ring  and  stringer  stiffened 

configuration.  In  addition,  the  geometries  were  chosen  in  such  a way 

that  some  comparison  can  be  given  with  stiffened  configurations  of  the 

only  stringer  or  ring  stiffened  type.  These  examples  correspond  to  a 

geometry  for  which  Z = 95.4,  (e  / 1)  = + 6,  (e  /t)  = + 3,  \ = 0.455, 

x — y xx 

p^_  = 100,  and  = 20  (see  Table  1).  The  imperfection  shape  is  taken  to 
be  similar  to  the  buckling  mode  and  the  imperfection  amplitude  is  varied 
from  half  to  four  times  the  thickness  of  the  skin.  The  reader  is  reminded 
that  the  analysis  is  based  on  the  smeared  technique.  According  to  Ref.  9 
this  configuration  without  rings  is  highly  sensitive  to  geometric  imper- 
fections when  the  stringers  are  positioned  on  the  outside  and  virtually  ] 

insensitive  for  inside  positioning  of  the  stringers.  The  present  results 
(see  Table  1)  indicate  that  with  the  presence  of  rings  the  sensitivity  is 
reduced  for  external  stiffeners  and  aggrevated  (increased)  for  internal 
stiffeners.  Examples  12-16  correspond  to  external  stiffening,  while 
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examples  17-21  to  internal  stiffening.  The  configuration  is  checked  for 
a primarily  axisymmetric  imperfection  (examples  12  and  17)  as  well  as 
symmetric  imperfections  (remaining  examples).  The  effect  of  the  imper- 
fection amplitude  is  checked  for  symmetric  imperfections  and  both  exter- 
nal and  internal  stiffening.  This  effect  is  shown  graphically  on  Fig.  3. 
Fig.  4 shows  the  plots  of  load  versus  "unit  and  shortening"  for  three 
amplitudes  of  the  symmetric  imperfection  and  external  stiffeners.  Table  2, 
which  corresponds  to  example  14  but  is  typical  for  all  examples  considered, 
depicts  in  tabular  form  the  computational  procedure  required  to  arrive  at 
the  final  results  shown  in  Table  1.  From  this  table  one  can  see  that  the 
minimum  total  potential  corresponds  to  n = 4.  For  this  n value,  all 
results  are  the  same  for  NP  = 35  and  NP  =65.  In  addition,  when  K =1 
and  K = 2,  the  critical  load  differs  by  4%  or  less.  Table  3 presents  a 
comparison  of  response  (radial  displacement  and  stress  function  amplitudes) 
and  critical  loads  for  different  K values  (number  of  Fourier  terms)  and 
two  different  values  of  NP  (number  of  grid  points). 


33 


She  1 1 . 


o 


r 


N 

XX 

^xx 

N*££ 

n 

k 

Points 

NP 

Total 

Potential 

"end 

shor tening" 

No . of 
Iterations 

CPU  TIME 
(SEC) 

19000 

0.54 

3 

1 

35 

-10777. 

0.021870 

3 

9 

22000 

0.62 

3 

1 

35 

-14499. 

0.025424 

4 

11 

23500 

0.67 

3 

1 

35 

-16590. 

0.027394 

4 

12 

25000 

0.  71 

3 

1 

35 

over  li 

mit  point 

13000 

0.37 

4. 

1 

35 

- 5028. 

0.015124 

3 

9 

16000 

0.45 

4 

1 

35 

- 7628. 

0.018620 

3 

9 

19000 

0.54 

4 

1 

35 

-10785. 

0.022161 

4 

11 

22000 

0.62 

4 

1 

35 

-14534. 

0.026189 

6 

20 

22750 

0.65 

4 

1 

35 

over  limit  point 

19000 

0.54 

4 

1 

65 

-10784. 

0.022160 

4 

19 

22000 

0.62 

4 

1 

65 

-14534. 

0.026188 

6 

33 

22750 

0.65 

4 

1 

65 

over  limit  point 

19000 

0.54 

4 

2 

65 

-10789. 

0.022232 

4 

70 

20500 

0.58 

4 

2 

65 

-12590. 

0.024197 

5 

90 

22000 

0.62 

4 

2 

65 

over  limit  point 

19000 

0.54 

5 

1 

35 

-10745. 

0.022248 

4 

11 

22000 

0.62 

5 

1 

35 

-14443. 

0.025847 

4 

11 

22  750 

0.65 

5 

1 

35 

-15459. 

0.026813 

5 

15 

23500 

0.67 

5 

1 

35 

over  limit  point 

22000 

0.62 

6 

1 

35 

-14362. 

0.025764 

4 

11 

22000 

0.62 

7 

1 

35 

-14318. 

0.025762 

3 

9 

22000 

0.62 

12 

1 

35 

-14267. 

0.025758 

2 

7 

35 

L 


r 


3b 


A 


limit  point 


0 


1 


2 


3 


4 


Fig.  2 


L 


Tx  8 


Radial  Dispalcemcnt 

along  y = 0 (Example  2) 


37 


Chapter  V 


CONCLUSIONS  BASED  ON  THE  NONLINEAR  BUCKLING  ANALYSIS 
A methodology  for  the  buckling  analysis  of  imperfect,  thin,  circular, 
cylindrical,  stiffened  shells  under  uniform  axial  compression  and  for 
various  end  conditions  is  presented  and  demonstrated  through  a number  of 
examples.  On  the  basis  of  the  results  reported  herein  very  few,  if  any, 
general  conclusions  can  be  drawn.  Among  these  one  may  list  the  following: 

(1)  It  seems  that  the  imperfection  sensitivity  of  generally  stiffened 
configurations  strongly  depends  on  the  curvature  parameter. 

(2)  From  the  examples  considered,  it  appears  that  configurations  with 
external  stiffening  are  more  imperfection  sensitive  than  configurations 
with  internal  stiffening. 

(3)  The  few  examples  considered  seem  to  support  the  contention  that  the 
most  severe  shape  of  imperfection  is  that  which  resembles  the  buck- 
ling mode. 

(4)  The  curve  that  demonstrates  the  effect  of  the  imperfection  amplitude 
on  the  critical  load  seems  to  be  approaching  a finite  asymptote 
(see  Fig.  3) . 

(5)  The  presence  of  both  stringers  and  rings  in  a configuration  alters 

the  conclusions  regarding  the  effects  of  positioning  of  the  stiffeners 
and  of  the  curvature  parameter  on  the  imperfection  sensitivity  of  a 
configuration  with  either  strings  or  rings  only.  The  limited  generated 
data  suggests  a nonlinear  coupling  of  the  individual  effects. 

(b)  In  general  one  should  not  expect  the  wave  number,  n,  corresponding  to 
the  limit  point  for  an  imperfect  cylinder,  to  be  the  same  as  the  one 
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predicted  by  the  linear  analysis  of  the  corresponding  perfect  geometry. 
This  difference,  though,  seems  to  be  minimized  when  the  configuration 
is  stiffened  in  both  directions  (see  Table  1;  Examples  12-21). 


IT 
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Chapter  VI 


IMPERFECTION  SENSITIVITY  OF  OPTIMAL  CYLINDERS 
The  formulation  of  minimum  weight  design  or  optimization  of 
stiffened  cylinders  under  destabilizing  loads  is  based  on  classical 
(linear)  buckling  analyses  (see  References  18  and  23  and  references  therein). 
A knockdown  factor  is  used  to  account  for  the  fact  that  such  geometries 
are  sensitive  to  geometric  imperfections.  Such  formulation  and  design 
procedure  is  fully  described  and  demonstrated  through  a number  of 
examples  in  References  18  and  23. 

The  precise  statement  of  the  problem  considered  in  the  above  two 
references  (in  its  most  general  form)  is  as  follows:  Given  an  internally 

stiffened  thin  circular,  cylindrical  shell  of  specified  material,  radius, 
and  length,  find  the  size,  shape,  and  spacings  of  the  stiffeners,  and  the 
thickness  of  the  skin  such  that  the  resulting  configuration  can  safely 
carry  a given  axial  compressive  load  with  minimum  weight.  Since  this 
configuration  is  sensitive  to  geometric  imperfections,  a true  solution 
can  be  accomplished  by  incorporating,  in  the  design  procedure,  a nonlinear 
buckling  analysis  (of  geometrically  imperfect  cylinders  - as  described  in 
Chapters  II  and  III).  This  task  is  formidable  and  it  will  probably 
require  an  unreasonably  large  amount  of  computer  time. 

The  reasonable  alternate  to  the  above  approach  is  to  follow  the 
design  procedure  employed  in  References  18  and  23,  establish  the  ODtimum 
point  in  the  design  space  and  then  perform  an  imperfection  sensitivity 
analysis  (nonlinear  buckling  analysis)  on  the  optimum  point  as  well  as 
the  surrounding  design  space  in  order  to  establish  what  effect  the 
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different  design  variables  have  on  the  knockdown  factor.  The  design  space 

that  needs  to  be  investigated,  for  a given  amplitude  of  the  imperfection 

(the  shape  can  either  be  considered  the  same  for  all  design  points,  or  it 

can  be  taken  to  be  as  one  that  leads  to  the  largest  reduction  in  load  - see 

Chapter  IV),  should  be  established  apriori  by  employing  the  following 

criterion:  if  for  a given  problem  the  optimum  weight  is  500  lbs.,  as  one 

moves  away  from  this  optimum  geometry  (variations  in  the  design  variables) 

the  weight  increases;  in  addition,  it  is  well  accepted  that  stiffened 

cylinders  are  not  as  sensitive  as  unstiffened  ones  (on  the  basis  of 

available  experimental  data);  therefore,  if  for  a given  imperfection  the 

critical  load  is  60%  of  the  corresponding  perfect  geometry  critical  load, 

(knockdown  factor  = .6),  then  the  design  space  in  which  comparison  studies 

are  made  should  be  such  that  its  boundary  weight  is  no  larger  than,  say 

50/.,  of  the  optimum  weight;  this  way  one  can  find  out  if  whatever  is 

gained  by  being  at  an  optimum  point  (on  the  basis  of  linear  buckling 

analysis)  is  not  lost  because  of  large  variations  in  the  knockdown  factor 

as  one  moves  away  from  the  optimum  design  point.  The  only  question  in 

this  procedure  is  whether  one  needs  to  use  one  value  of  the  load  and 

one  knockdown  factor  in  the  linear  optimization  procedure  or  several  in 

order  to  establish  how  the  optimum  design  point  moves  In  the  design  space 

with  small  variations  (10-15%)  in  the  N value  used  in  References  18 

xx 

and  23.  This  question  is  dealth  with  in  the  section  entitled  "Recommended 
Design  Procedure." 


43 


1 . Design  Examples 

This  alternate  approach  is  employed  herein  and  it  is  applied  to  the 
following  two  design  cases  of  References  18  and  23. 


Case  1.  R = 95.5  in.,  L = 291  in.,  N = 800  lb/in., 
xx 


(Ref.  18) 


E = 10.5  x 106  psi;  p = 0.101  lbs/in'?,  v = 0.33, 


6 = 0.0442  in.  (imperfection  amplitude) 

RSRR  (rectangular  stringers  and  rings) 

MG  = 0.02  in.  (minimum  gage  constraint) 

Boundary  Conditions  SS3. 

Case  2.  R = 85  in.,  L = 100  in.,  N = 2700  lbs/in., 

xx 

(Ref.  23)  E = 1Q<5  psi>  p = 0101  i bs/ iru  , v = 0.33 

6 = 0.1  in.,  RSRR  (rectangular  stringer  and  rings), 

TSRR  (T  stringers  and  rectangular  rings  with 

C = 1.079  and  1.135) 
x 

MG  = 0.05  in..  Boundary  Conditions  SS3. 

in  both  cases  the  imperfection  amplitude  is  almost  twice  the  skin 

thickness  of  the  optimum  geometry.  Both  cases  are  checked  for  an  axi- 

symmetric  and  for  a symmetric  shape  of  the  imperfection.  The  results  for 

the  latter  shape  are  reported  herein,  because  they  yield  the  greatest 

sensitivity.  This  shape  is  given  by 

o . . . mnx  ny 

w (x,y)  = 6 sin  — — cos 

Case  1 corresponds  to  a geometry  with  rectangular  stringers  and  rings 
(RSRR)  and  the  optimum  design  geometry  is  reported  in  References  18. 


Case  2 corresponds  to  a geometry  with  tee  stringers  and  rectangular  rings 


L 


(optimum  stiffener  shapes)  with  various  flange  width  to  web  height  ratios 
(k  = 0.65,  C = 1.212;  k = 0.45,  C = 1.135;  k = 0.30,  C = 1.079) 

X X X X X X 

including  the  limiting  case  of  k^  = 0 and  C^_  = 1.000  which  corresponds 
to  rectangular  stringers.  The  corresponding  geometries  are  reported  in 
Table  6 of  Reference  23. 

The  results  are  presented  in  tabular  form  and  are  discussed  below. 

2 . Discussion  of  Results 

Case  1 . The  important  results  associated  with  the  investigation  of 

Case  1 are  presented  in  Table  4.  The  minimum  weight  configuration  obtained 

in  Reference  18  is  labeled  as  point  2 on  this  table.  The  knockdown  factor 

for  this  point  is  found  to  be  0.9125  and  the  shape  of  the  imperfection 

that  yields  the  greatest  reduction  corresponds  to  the  perfect  geometry 

buckling  mode,  i.e„,  m = 17  and  n = 9.  Because  of  this,  the  imperfection 

for  all  design  points  in  the  neighborhood  of  the  optimum  is  taken  to  be 

0.0442  sin  12l2L  cos  and  n for  each  point  is  a free  parameter  and 
L R 

evaluated  such  that  it  yields  the  most  reduction  from  the  classical 
buckling  load.  Points  1,  3,  and  4 correspond  to  the  optimal  geometries 
for  spanning  values  of  the  curvature  parameter  Z . Note  that  as  Z decreases 
the  knockdown  factor  decreases,  which  means  that  the  thicker  the  shell 
the  bigger  the  imperfection  sensitivity.  By  comparing  the  weights  and 
the  knockdown  factors  for  these  four  points,  it  can  easily  be  concluded' 
that  point  2 geometry  is  better  than  the  geometry  corresponding  to 
points  3 and  4.  The  only  question  exists  with  point  1 geometry  since  it 
is  less  sensitive  than  that  of  point  2.  But,  since  the  weight  of  point 
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1 geometry  is  heavier  than  that  of  point  2 by  5.5%  while  the  knockdown 
factor  is  only  higher  by  2.8%,  one  might  still  conclude  that  the  optimal 
geometry  (point  2)  is  the  best  even  in  the  presence  of  imperfection  sensi- 
tivity. At  this  point  the  authors  would  like  to  point  out  that  these 
conclusions  do  not  suffer  by  keeping  m = 17  in  the  imperfection  shape, 
because  points  3 and  4,  only,  might  be  affected  by  letting  m be  a free 
parameter,  in  which  case  the  knockdown  factor  for  these  points  might 
become  smaller. 

Once  the  effect  of  the  curvature  parameter  is  established,  the  inves- 
tigators considered  the  effect  of  the  other  design  parameters,  such  as 


extensional  and  bending  stiffnesses  (a  ,a  ,\  ,A  )•  For  the  Z value 

° x y xx  yy 

that  corresponds  to  the  optimal  geometry  (point  2)  the  design  space 


surrounding  the  optimal  point  is  investigated  (in  the  presence  of  imper- 


fections) and  the  comparison  of  the  weights  and  knockdown  factors  (points 


2a,  2b,  2c,  2d,  2e  and  2)  supports  the  conclusion  that  point  2 yields  the 


minimum  weight  design  even  when  imperfections  are  present.  In  this  com- 
parison it  is  also  observed  that  small  variations  in  all  other  design 
parameters,  except  the  curvature  parameter  Z,  have  a very  small  effect  on 
the  knockdown  factor.  A similar  investigation  of  the  surrounding  design 
space  of  point  1 (Z  = 41,850)  yields  the  same  conclusion  (these  results  are 
not  reported  herein). 


Case  2.  The  results  of  the  investigation  for  this  case  are  reported  in 

Table  5.  The  imperfections  are  takens  as  A)  w(  = 0.1  sin  cos  and 

L R 


The  optimal  geometry  corresponds  to  tee 


stringers  and  rectangular  rings  and  is  labeled  as  point  1 in  Table  5 (see 
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Table  6 of  Ref.  23).  Two  shapes  for  the  imperfection  are  considered  , 
one  denoted  by  0.1  sin  — cos  and  the  other  by  0.1  sin  ^7-^  cos  — 
where  m in  the  second  is  taken  to  be  the  classical  buckling  mode  number 
of  half  sine  waves  in  the  axial  direction  and  n (in  both  cases)  corres- 
ponds to  the  number  of  full  waves  in  the  circumferential  direction  that 
yields  the  smallest  buckling  load  (in  the  presence  of  the  imperfection). 
The  imperfection  amplitude  in  both  cases  is  twice  the  thickness  of  the 
optimum  geometry.  Point  2 corresponds  to  the  optimum  geometry  at  a 
lower  Z value.  For  both  points  (1  and  2)  the  geometry  corresponds  to  tec 
stringers  and  rectangular  rings  (TSRR)  . Point  3 employs  rectangular 
stiffeners  and  corresponds  to  the  optimum  geometry  for  this  construction 
and  the  same  Z as  point  1.  Points  la  through  lg  correspond  to  the  same 
construction  (TSRR)  and  Z value  as  point  1 but  different  values  of  the 
remaining  design  parameters  (_a^,a  ’^xx’^yy * Finally,  points  3a 
through  3c  correspond  to  the  same  construction  (RSRR)  and  Z value  as 
point  3 but  different  values  of  the  remaining  design  parameters. 

Optimal  geometries  corresponding  to  higher  Z values  than  point  ; 
are  checked,  but  since  the  knockdown  factor  for  these  geometries  is 
smaller  thar  that  of  point  1,  the  results  are  not  presented  herein.  This 
conclusion  is  the  same  as  that  of  Case  1. 

In  making  the  comparison  studies  the  smallest  of  the  knockdown 
factors  is  considered. 

In  comparing  points  1 and  2,  the  results  are  inconclusive  because 
the  difference  in  knockdown  factors  and  that  of  the  corresponding  weights 
arc  approximately  the  same.  This  might  be  interpreted  as  both  geometries 
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being  very  close  to  the  optimum  for  this  construction. 

A comparison  of  the  data  for  points  la-lg  with  those  of  point  1 
suggests  that  point  If  might  yield  a lighter  conf igurat ion  than  that 
of  point  1.  The  reason  for  this  conclusion  is  based  on  the  observation 
that  the  geometry  of  point  If  is  heavier  than  that  of  point  1 by  IX  while 
it  is  less  sensitive  by  approximately  15X. 

A similar  comparison  among  the  RSRR  geometries  shows  that  point  3 
geometry  is  indeed  the  best  geometry  (same  conclusion  as  in  Case  1). 

Finally,  a comparison  between  the  two  different  constructions 
suggests  that  the  best  of  either  corresponds  to  an  acceptable  design. 

The  only  thing  in  favor  of  the  RSRR  construction  is  that,  if  cost  of 
manufacturing  is  taken  under  consideration  (additional  constraint), 
then  this  construction  is  superior. 

Note  that  in  both  cases  the  results  of  the  linear  optimization 
procedure,  corresponding  to  only  one  value  of  N , are  employed.  This 
may  cast  some  doubt  to  the  validity  of  the  conclusion  and  for  this  reason 
an  improved  design  methodology  is  suggested  in  the  next  chapter. 


48 


Table  5 . Effect  of  Imperfections  on  the  Optimal  Ceornetry  (CASE  2). 

o - . 5ttx  ny  n o . mrrx  ny  1 

A;  w 0.1  sin  cos  ~r~  B;  w 0.1  sin  cos  I 

_ L L Li  R 
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0.05 
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35 
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0.5386 

0.3310 

0.4608 

0.1351 

0.2128 

0.2309 
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40.04 

66.35 

165.45 

430.88 

831.17 
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0.354 

0.325 
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7 
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0.8541 
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0.05 

15 

40 
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0.05 
15 
35 

0.4940 

0.1804 

111.14 

220.95 
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484 

0.05 

17 

35 

0.5990 
0.1105 
173.11 
135.38 
0.483 
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1.079 
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0.05 

13 

40 

0.4525 
0.2685 
76.47 
429.66 
0.3  76 
1.025 

1.079 
2221 
2 700 


Chapter  VII 


RECOMMENDED  DESIGN  PROCEDURE  FOR  AXIALLY  LOADED 
IMPERFECT  STIFFENED  CYLINDERS 

As  is  seen  from  the  discussion  of  the  results  of  the  two  design 
cases  considered,  there  are  no  general  trends  and  conclusions  that  can  be 
applied  to  all  design  problems.  In  addition,  the  entire  study  (presented 
herein)  and  the  resulting  conclusions,  as  far  as  the  position  of  the 
optimum  design  (including  the  effect  of  geometric  imperfection)  in  the 
design  space  is  concerned,  are  based  on  the  following  conjecture:  as  the 

applied  load  N (including  a knockdown  factor  - used  in  the  linear  opti- 
mization procedure  of  References  18  and  23)  experiences  small  changes 
the  corresponding  linear  optimum  point  shifts  smoothly  and  by  small 
amounts  in  the  design  space.  For  example,  if  for  a given  problem  the 


value  of  N changes  from  800  to  880  (107.  increase)  the  corresponding 
xx 

optimum  design  is  expected  to  change  from  Z = 35,000,  a ~ 20,  q>  = 90, 

x y 

X = 0.64  and  X = 0.20  to  Z = 35,000  + AZ , cr  = 20  + Aa  > etc.  where 
xx  yy  x x 

the  deltas  denote  changes  of  less  than  107,  from  the  previous  values.  To 
further  clarify  this  point  and  its  implications  consider  the  Case  1 design. 


First  of  all,  the  optimum  (linear  optimization)  design  is  characterized 


by  (point  2 of  Table  4)  h = 0.0221  in.,  q-  - 20,  a ~ 95,  X - 0.6319 
J r x y xx 

and  X = 0.2049  with  a weight  equal  to  755  lbs.  The  value  of  the  applied 

yy 

load  used  is  800  lbs/in.  Since,  at  this  geometry,  the  knockdown  factor 


is  0.9125,  then  the  shell  is  in  reality,  designed  to  carry  safely,  with 

minimum  weight,  a load  of  730  lbs/in.  (N  = [730/knockdown  factor]  = 800). 

xx 
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Realizing  that  the  knockdown  factor  is  not  known  apriori  and  observing 

that  this  factor  depends  on  the  particular  point  in  the  design  space, 

then  different  values  of  N must  be  considered  in  the  entire  optimization 

xx 

procedure,  before  one  is  able  to  zero  in  to  the  final  optimum  design  point. 

Because  of  this,  the  following  procedure  is  suggested  in  order  to 
achieve  the  minimum  weight  design  for  a given  applied  load  (for  the  sake 
of  argument,  say  1000  Ibs/in.)  in  the  presence  of  geometric  imperfections. 

a)  Specify  the  maximum  allowable  amplitude  of  the  geometric  imper- 
fection. 

b)  Guess  the  value  for  the  knockdown  factor  (say  0.8  for  the  ensuing 
discussion)  . 

c)  Employ  the  design  procedure  of  References  18  and  23  with  N = 

(1000/0.8  = 1250,  and  find  the  minimum  weight  design  point 

(W , o , a , 1 and  A ) . 
x y xx  yy 

d)  If  the  actual  knockdown  factor  is  within  10%  from  the  guessed  value 
of  0.8,  then  perform  comparison  studies  as  outlined  herein. 

e)  Repeat  steps  c)  and  d)  for  values  of  N equal  to  1000/(0.72) 
and  1000/(0.88).  An  examination  of  the  three  comparison  studies 
will  yield  the  optimum  configuration  within  any  desired  (reasonable) 
accuracy. 

f)  If  tlie  actual  knockdown  factor  differs  by  more  than  10%  from  the 
guessed  value  (say  20%,  lower),  then  repeat  step  c)  for  values  of 

N equal  to  100/(0.55),  1000/(0.65),  and  1000/(0.80)  and  perform 

xx 

the  comparison  studies  as  mentioned  in  step  e) . 

By  following  the  above  procedure,  the  designer  is  assured  of  arriving  at  a 


53 


minimum  weight  design,  which  can  carry  safely  a uniform  axial  compression 


Appendix  A 


BUCKLING  OF  IMPERFECT  STIFFENED  CYLINDERS  UNDER  COMBINED 
AXIAL  COMPRESSION  AND  PRESSURE 

The  buckling  analysis  of  geometrically  imperfect,  stiffened,  thin, 
circular,  cylindrical  shells  under  uniform  axial  compression  has  been 
presented  in  Chapters  II- IV.  The  purpose  of  this  Appendix  is  to  show  the 
extension  of  the  analysis  for  the  load  cases  of  hydrostatic  pressure, 
lateral  pressure,  and  combined  uniform  axial  compression  with  lateral 
pressure . 

Most  of  the  investigations,  reported  in  the  open  literature  which 

deal  with  buckling  of  imperfect  configurations,  are  for  uniform  axial 

compression  and  isotropic,  constant  thickness  geometries.  Hutchinson 

9 

and  Amazigo  investigated  the  buckling  of  imperfect  stiffened  cylinders 
(stringer  and  ring  stiffened  only)  under  hydrostatic  pressure.  There  is 
no  reported  investigation  for  stringer  and  ring  stiffened  geometries. 

The  self-equilibrated  lateral  load  p(x,y)  [positive  inward]  is 
expanded  in  terms  of  Fourier  series 

K 

P(x,y)  = ^ p^x)  cos  (A  - 1) 

i=0 

Then  all  equations  in  the  body  of  the  report  must  be  changed  appropriately 
to  reflect  the  inclusion  of  this  load.  For  example,  the  term  -po(x)  must 
be  added  to  F.q . (26a),  the  term  -p^(x)  to  Eq.  (26b),  and  the  following 
term  must  be  added  to  the  total  potential  expression,  Eq.  (27). 
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K 

-ttR  ^ 2p  (x)  W (x)  + p.(x)  W.(x)  dx 

i_  rov  ' o U x l j 

i=l 


The  related  computer  program  (Appendix  B)  includes  these  changes. 

A number  of  illustrative  examples  are  considered  and  the  results  are 
presented  in  tabular  form  (see  Table  A-l)  for  uniform  pressure  only.  For 
all  examples,  the  boundary  conditions  are  taken  to  be  the  classical  simply 
supported  ones,  SS3.  The  program,  though,  is  written  for  all  types  of 
lateral  (clamped,  simply  supported,  free,  etc.)  as  well  as  in-plane 
boundary  conditions.  Some  of  the  examples  considered  serve  as  bench  marks 
for  the  developed  methodology  and  computer  program.  In  addition,  new 
results  are  reported.  For  all  cases,  two  types  of  geometric  imperfection 
are  considered  which  are  characterized  by  § = 1 and  | = 2.  These  are: 


5 = 1 


o . . mrrx 

w (x,y)  = t sin  -j— 


+ O.lt 


s in 


cos 


tV£ 

R 


(A -2) 


p = ? 


w°(x,y) 


mnx  ny 

t sin  — — cos  -r- 
L K 


(A-3) 


Note  that  the  imperfection  characterized  by  Eq.  (A-2)  virtually  denotes  an 
axisymmetric  type  of  imperfection,  while  the  one  characterized  by  Eq.  (A-3) 
denotes  a symmetric  type  of  imperfection.  In  both  cases  the  amplitude  is 
approximately  one  skin  thickness.  The  values  of  m and  n that  correspond 
to  the  most  sensitive  geometry  are  reported  in  Table  A-l. 

The  examples  reported  have  the  following  common  geometry 


L = 4";  R = 4";  t = 0.04" 
v = 0.3;  E = 10.5  X 106  psi;  Z = 95.3 


(A -4) 
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Note  that  examples  1 and  2 correspond  to  internal  ring  stiffening  only, 
examples  3 and  4 to  internal  stringer  stiffening  only,  examples  5,  6, 

8-11  to  internal  ring  and  stringer  stiffening,  and  example  7 to  external 
ring  and  stringer  stiffening.  In  addition,  the  loading  for  all  the  examples 
is  lateral  pressure  for  examples  1,  3,  and  5,  hydrostatic  pressure  for 
examples  2,  4,  6,  and  7,  axial  compression  for  example  8,  and  combined 
loading  for  examples  9-11  (internal  pressure  of  half  and  one  atmosphere 
respectively  for  examples  9 and  10,  and  external  pressure  of  one  atmos- 
phere for  example  11). 

Examples  2 and  4 are  taken  from  Ref.  9 and  the  results  are  in  very 
good  agreement. 

Among  the  important  observations,  on  the  basis  of  the  generated  data, 
one  may  list  the  following: 

1)  For  pressure  loading  the  ring  only  stiffened  cylinder  is  sensitive 
to  geometric  imperfections  while  the  stringer  only  stiffened 
cylinder  is  not.  This  observation  is  true  for  lateral  pressure 

(N  = 0;  examples  1 and  3)  as  well  as  hydrostatic  pressure 
(N  = pR/2;  examples  2 and  4).  Furthermore,  the  presence  of  the 
axial  load  in  the  case  of  hydrostatic  pressure  (compare  results  of 
examples  1 to  2 and  3 to  4)  aggravates  the  imperfection  sensitivity. 

2)  When  stiffening  in  both  directions  is  used  the  configuration  is 
still  sensitive  but  not  as  sensitive  as  the  ring  only  configuration 
(compare  the  results  of  examples  5 to  1 and  6 to  2). 

3)  It  is  concluded,  in  Ref.  9,  that  positioning  the  rings  on  the 
outside  (rings  only  configuration)  makes  the  shell  more  sensitive 
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to  geometric  imperfections  as  compared  to  internal  positioning 
of  the  rings.  This  observation  is  also  true  for  stringer  and 
ring  stiffened  cylinders  (compare  the  results  of  example  7 to  6) . 

4)  For  stiffened  cylinders  (both  rings  and  stringers)  under  uniform 
axial  compression  with  or  without  lateral  pressure  the  presence 
of  internal  or  external  pressure  of  up  to  one  atmosphere  does  not 
have  an  appreciable  effect  on  the  final  results.  Linear  theory 
predicts  minimal  effect,  in  the  right  direction  (examples  8-11) 
(negative  p implies  internal  pressure)  and  nonlinear  theory  shows 

that,  for  the  same  imperfection,  the  sensitivity  is  not  affected 
by  the  presence  of  the  small  internal  or  external  pressure. 

5)  The  buckled  shape  in  the  circumferential  direction  (n)  is  vir- 
tually the  same  for  the  nonlinear  theory  in  the  presence  of 
geometric  imperfections,  as  it  is  for  the  perfect  geometry  linear 
theory . 

Finally,  note  that  all  results  are  derived  for  one  value  of  the  cur- 
vature parameter  Z and  two  stiffener  geometries;  therefore,  no  generality 
should  be  attached  to  the  above  observations.  If  one  is  interested  in 
performing  extensive  parametric  studies  in  order  to  draw  general  conclu- 
sions the  present  computer  program  is  very  efficient  (Appendix  B) . 
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Table  A-l.  Buckling  Results  of  Imperfect  Stiffened  Cylinders 
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COMPUTER  PROGRAM 
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Z)  COhtoON/BOUND/LM  ,LSH 

W> ovAVvdlOo^  CowA'\U»v\  oA  4tvi?  (^LS  i.), 

ocV.  <Ua  ^\v^CtS^  o\  4U  SWW. 

%)  Common  /FI  D F ft  / bE  tTA  , ALl  , & ft  1,  A V.2  jfeT2,Cft2  • 

Coe^ioevts  o^-  ^Av\\te  tiv^e^e^oe  |«rw\  , ^ 

f Z V^A,  «^=^rl/ZSf  >f*‘V£  • 

4)  Co(v\tAOM/FQ^gxe./  vc^obR  , ^ ^4,^3,  V£,  ■ 

Fovx^e^o  server  ^v*VL  (vc  = VC*mO  ?o^w^ters 

c^£^^)€/v\^£a^A-  O vv  ft.  , 

5)  I / £ & ,w , VU-i 

6)  Cot^ON  /rA,CTOe./C-L,CZ,  - - - -C12  ■ 

Coe^'c*e*vU  uiftvd  <xv«  U^amUv^  ca  Circus  \&re^U  ujo^/e  *«~W. 

7)  Co^rwON/rACT2/t)\.i-p\_4,Dftl-DA4;OSZ>,PB>3>p^A>yNI>  1 

Coe|^»c\eAAts  ^ ^ /3-j  ^4  > ^°i  > W , \r>  Exxp  . 
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8)  Common  /cpisk/x2i  (soi)/i  Csoi  ) 

Dtvtcl-  access  dcjto^  set  Zi  o^l 

Cofv>moN/PACT  3 /PIS  , XL ,XH  . 

^&*^vv\ei£'r/  XL-L  , X H = -fc  . 

V0)  Coiv\iwoM/pReS  t / VVN\  (2Q0;5)^ETfV\  fgQOy5Xw(WPUO0>') 
TVv€  vecLo'c  o^.  )preoio^s  so^wtua-vv  , a.wA 

aA-  ^ov  fo^veA-  "tje/cws  C . 

lij  CoN\moN/PRes2/w^(aoo,s),WZP(zoo;5^w,z.PP  (•zoo  ys) 
Iws^^ecl^vN  v^c<Le  W°  , W0'  W°"  oX  po'iwl  1 
'^I'T  Aooorve^*  "tiircvvv  L , 


*2)  Common /pres  3/  Ffy\  (zoOjftj^yf  rv\  (^oOjB)  ;^(v\P(eoc;8) 
A(~e.  uedox-  p-reoiovAS  so5LoVvxv%  ^ ^ ,avN<J.  ^ 

&A  f>CUv\t  1_  ^<nr  'Tov^.'rie^'  te*-™  t . 

13)  COMMON  / XxU?At)/XPR£-S 

X?R t S - k^ArosteA'c  ^?r€SsiAve  . 
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2. Flow  Cvia^t  for.  pr6pp>>p.a.t\OH 


STAEf 


/title 

9A8 

Av\y  al|pUckvw.vv\«.Y\c  y vx^orwiOvlUjovN  S^K 
Ox.se  (_ 72  oo  5_lcw^  s ^ • 

/ N£QPOT  VCFOUR  ,LSlj  LSN,V/RlNT  ,LIY\OD, 


*PP£T. 

N 6QP0^  ;No.  o(  ?o\*ts  \ * X-lir«<^urv\  ^ (^s<  NtQ?OT^"EOoy 
\OF01AR  • ^ob(i(tR  S£E>ts  i-sOMT  (r')  ) ^ O-UFCVift  <4^). 
It  is  POSSIBLE  To  INCREASE-  KPOUE  Mo  N6Gl?0T  fey 
CUANCrlN'C*.  D\W\E(sfSvON  KNO  COW^ON  STfWcfAtNT  . 

Vs ^ ' Boundary  cdmP\T\ok)  op  the  F'P-st  Po»vjt. 
USN  '.  BoUNOA«.Y  CONDmON  OF  "THG  L-ftST  PoiWTi 
Note-  ».  LS  = 1.  { ov  t L^-2  ^ SSZ  ) LS  = 3>  ^ 

LS  -4  fcr  SS4-  ) LS--S  4^  CCt  , V-S-G  ^ CC^ 
LSr7  -fiv  CC3  , LS-g  ^«v  CC4,US-^  FFA-; 

\LS'10  |V  5'jvrtinetirvj  Cow^vtce-vi  LS  - i-A-  '|*v 

s>Avv\*vvet^'-A  Coa\A.\4,\ova  . \ 

^ ^ (Gow-tU) 


(Ccmt-’d.) 

0 - w'>v\»v«(aa»v  pr'V\\-oict  < ^ r t*».x\yvha.vw  yr > v&e vxt 

LOAOO  0 - <1°^  mi  privd  vnoAto  , * : pv'w^  wxO&es 

X 0 D tT  0 Ao(r>  v\*A.  UX  C <1  eij^W'vV' v v\ OiAxt 

^ 4eiLeA'vv\>v\oUv^-  ^)r\v\fc-i  it  . 

f\?  (l2Wi4-,n^4)  Cp(ttvc+4.,vafc4  4'i 

?£  (l2W+A;iX^-v4)  , GP  (UM4;/L)>X?(l2V:+4;l^TX.(nWX4\C.02i!:v4)y 
W\T  (\2K+4)  ,Vl(0*k-v4)*(>2K44)). 
cowwcn/ps^V  ww  (weft?oT^O,  exf\  (wea?<rr  ,e*0,yvn\P  (h  tcaPox;  MO  . 
CO|y\MON/pR£S  2/  vv/z  (NEQPCiTjM  l}  NJZ?  (NtGtf  (H  ,V^;>NZ?P(NtQ  ?CT\  i) . 
C0<v\M0V/pfcts3/p  tv\  (wE(XPOT^vc)^F(y\  (*iefcPcx>0,w(NEG?ox  /zk). 


Tk«  lA5«r  mmIl  svJk'co\Aw't\Y\e.  T(V\PtRP  ^-<jy -VU*. 

4Us.  kW\^tvl^-e.ctiov\  <kmA  A&av7&\,’w/e  s . 


WZ  (x,T^  r VV°  (.  iposi  b\’ue  \ vxojck.'rj- ) 
WZ.pCljI)  -W° 
r \N°  " 

X --  1,  NE.QPOT  w\<?$U  ^powv^ 

J - 1 Vc^OUR-Vl  4°^  '^o^'er  'fce*-w\<> 


Mo\^-  ~HaX  Svvvce  yg0  = YJ^^Otos  — ^ 

v.  - o 

*-VC«Av\ 

T-l  4or  i'=°  , *~X  3^KV0UR  + 1 V tpouft 
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6 E 12.4 


DUJ^'l  ’For  ^ixe<t  (lejDml  ^Vis^rt  \ ^wl  NxCf 
Dt-KiC)=2.  for  ^■xei-  (Kx\oJL  covA'Qre^svow  V^xx  > \iv^  Vcr ' 
{)LNOr3  I^Xvoi  5~t>fcvA  \pcess\Ave  O.'ce 

V>^  Uu.  |^tW  XtfXX  ( VI*  -r*NX**'X?fctO 

XMXX  Toy  OLWQ  :1  x\  -\W  Xv^vajL  o^;»jL  do<A 

'Fcty  0^-N^>  - *2.  W 4&J.  ^ixe<S.  OjxvJI  ^otM^  • 

VtK  0 t KJ  O - 1>  Vi  'TW  ^JLrr  ALaX  'reiUiti 

N**  9vwX  ^ovUveXNKX  .-Cfiv^^aievi). 

xpaec,  Voy  OumO-1  va  4Ui  |ixe<l  ^Y«ss^ve 

^cn-  iXNT}  ~ 2.^  \<>  -V-L-ft.  WdtvoJl  OY'CViu.'C^ 


DNX9 


&CCUB. 


fcl 


( 'Pos » V-iue 


V V\V^J 


9oy  tJvWQ-i  ^o'vU  \v>tYew\^Ask.  m cwv«JL  (^o^ 

\-~oy-  v)u^Dr2)^  v»  vVi/v£  VV\Crt»V^^-  W>  ^>Y^S^IaY^ 
re^uvYeJl  jLo«jd  C^CCtv>rOvCv|  \v\  ^7  ^/x~C.£^v^- 

(V\  Ova  iwww  v\ vx.yv\Q)6/y  o^.  ioaA  >po\v\X^ 

* - XMXX  -V  9.1=#  ONX\>  (^yT)\.n0-:O 

y - R.I*DNX?  ^avOLNO-V). 


© 


1016 


TUe  iv\\k&d  Sokvtuw  -^«v 
'AowQ,vv\ecvr  ^^£/w\  wme-o 
XU.  C^wt  9*c  So^AaIa^vx 
it /L  XV\  vtuiA  ^CjIaAaXVx  ^v  ttoi 
v\ov\t\ v\t«.v-  ‘jvjsteA^  (^xWxx  A- t>wxP 
0Y  XVft-ti-t  T>^x^  ^)  ts  4iv«  SoLI^cm 
®dt  4U  ^r^uvcnM  SLooA  (ce. 

XNXX  0Y  X?£csV 

(^Ti  ^ Xc  -Uk*  iuow-1.^ 


XN  oh 


INON  - 1 


tNON  - 2- 
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5 


N1NN  } UNtSVI  ,ILNW 


1016 


fMNN  ~ UI<XW€.  ysva/wJ^y,  Y\  ,~rive  yroC^rtw 

^aaAs  4W  il*yv\i-)c  ipOwvfc  ^ov  "iXvvo  . 

LNNN-0  : ~TW  'proc^rcw^  doeo  ujWJ^  v%  (^oe<»  AWc. 

VY\  \ Y\  \ vwsvA-Wv  ^o^€m)uA^  eA\err^^ 

LMNN  -*i  Tk*  flrOipAvn  |^vv<U  AW  hvvuk  ^wfc  ^ ^NN, 
ca$cA.(todte.s  -AW  totoJjL  ^ot&v\UA£ 

SvAYOuv^.'\w.^  KiKtsj  ©A.  Qxjoc^L  !lav>e^'> 
§Q\>yw  iUwsv^  ^owvA-  (kjnn"). 

LNW  N ~ 2 TW  pro^vx>wvw  OX^CAj^ltx^o  AW  toWS-  ^C^AvOUaS 

CkQ/^  <^CfY  '4Wt.  (LoOw^  XKiXX  ."Tw  At\»Vo 

Oise  -Set  XNXX  ciloS€  W>AWl  Q/>tv  wvoWst-  Hiwwt  ^JO'vfc 

"l€v\  H '-W-Q.  YaUj>.  v wwava  VSwvw^-Q^/v  o V\  ~\i0^r 

{!aa£o  ujUacIv  A-We.  iprO^OvAAA  (^dLuJ^iw 

tcR/iJl  ^oW-Loil  v'vs  cprduuv  A*  ^AW«L 

WyA\vv/v*»~  [oa!  u>'rv'es^ov\(Lv<^  vn  V 
Vtfv  e/*  owA^9-e  , \ v\  = s <*/vJ)l  Ivn  'M  r.  'Z.  +W 

yn^vxx av\  UX.^.Caa&oAao  AW  \x)W5  yot^AjWil  1*^  fc-  ^ < C>. 

UT  C 0t  (v'-f ) -vWaa  i-V.  Ca^aJU^, 

LJt  v\-7  ,‘v|  vsot  ^Y  w-a, 


X Ul\N  ^ 10 
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lcontn 
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III. 


COMPUTER  PROGRAM 


72 


i 

I 


FRGGRAF  N AIMINPLT  , OUTPUT,  IAFF5=  INPUT,  TAPE6  = 0UTPUT,  TAPE2I  ,T  APE21  ^ 

1 1 T A P _ 2 c ) 

C POST  COCKLING  OF  SI1FFEMED  CYLINDRICAL  SHELLS  UNL^R  UNIFORM  AXIAL 
C COMPRESSION  (NONLINEAR  THEORY) 

COMMCN/CINTG/NEQPOT.MI  (5.0 
CCMtiCN/BCUNC/LSl  »L  SN 

CCMMCN/F10FR/DEL TA,AL1,GA1,AL2,BT2,GA2 
C0MMCN/F0URIR/KFCUP,K6,K4,K3,K2,K1 

CCMMCN/GE0M/RR,0C,Hll,H12,H22,Ull,Qi2,C22,Ull,C12,D22 
COMMoN/FACTCR/Cl ,C2,C3,C4,c5,C6,C7,Ca,C9,C10,CU,Cl2 
CCMMCN/F ACT2/DL1 , D L2 , 0 L3 , DL4 , C A1 , C A2 , D A 3 , 0 A A , 0 B 2 , DC  3 , 08 4 , XNl , E X XP 
CCMMON/COISK/I21 C5'l) , 122  <501> 

C0MM0N/FACT3/DL5  , XL, XH 

COMMON /PkE  SI /WM  < 2l‘  O 5 ) , ET M ( 2 0 f.  , 5 ) , WMP  ( 20  * , 5 ) * * 

C0MH0N/PRES2/WZ ( 20C , 5)  ,WZP  (cl C,5 ) ,WZPP  (2j  j ,5)  »» 

CONMUN/PEESJ/FK ( 2C  ~,8)  ,XFM  (20r ,6 ) ,FMP(2nC,8>  ** 

C GM.ION/X XL 0 A 0/  X P RE S 

DIMENSION  W W M ( 5 ) » F F M ( 8 ) ** 

DIMENSION  TI(10) 

DIMENSION  WF  (2,5)  , XWF  ( 2 ,5)  ,FF (2, 8 ) ,XFF  (2, 8)  ** 

DIMENSION  AF (52, 52) ,cp  (52,52) ,C° (52,52) , PR  ( 52 , 52 ) ,GP(52,i)  ** 

DIMENSION  XF  (52 , 1 ) , T1 ( E2>  ,CC  (5c)  , MT ( 52)  , VI (273  A ) * + 

C Z7lk-bc*b2 

DIMENSION  WCCN ( 2 0 , 5)  , FCON  (2C  ,8)  ** 

C ALL  THE  CARLS  WITH  SIGN  **  IN  COLUMNS  73,74  DEPENu  ON  NUMBER 
C OF  FCINTS  AND  KFODR 

EQUIVALENCE  ( A P ( 1 » 1 ) , V 1 (1 ) ) 

CALL  CFENMS (21, I 2i  ,5^ 1 ,n> 

CALL  CFENMS(22, 122  ,5C1  ,0) 

ECCNV= r . 0 uE CO 

Mi  A X N = 5 2 ** 

M A X 2 = F A X N * rl  A XN 
NR  H S= 1 

N J = 2 C n ** 

NW  = 5 ** 

NF  = 8 ** 

C NJ,Nw,NF  - FCR  DIMENSION  --  NJ-MAXlMuM  POINTS  IN  AXIAL  DIRECTION 
C NW=  MAX1MLM  KFOUR+1  , NF  = MAXIMUM  2*KF0UR  , M A Xi.  = lc*  KF  OUR+  4 
C IN  CkQ  R K INCREASE  THE  CAPABILITY  OF  THE  PROGRAM  FOR  MANY  POINTS  IN  AXIAL 
c OIRECIICN  AND  HIGHER  LIMIT  OF  F CUR  I ER  ScRIeo  ThL  USER  H aS  TO  CHANGE 
L ALL  THE  CA,<US  WITH  ThE  SIGN  **  IN  COLUMNS  7 3 AND  / 4 
1111  WRITE  (6,20 

RE6D(5,1L ) (TI(I)  ,1=1,5) 

WR  ITE  (E,  60 

Wk I T f (fc,lC)(TI(I),I=l,p) 

RE Au (5 ,10  O NEUPO  T,  KFOUK  ,LS1  , LSN  , i_FRI  NT  ,LM0U,1DDET 

IF  (LPF.INT.cC.l)LMCC  = l 

READ(5,2uORR,XL,XH,£LAS,XN1 

HEAD (5 ,2  DC ) XLAMD , YLAMC  , E X X , E » Y , R HOX , RHCY 

cX=-E>X 

l Y = - E Y V 


oALL  CDEFF(EX,EY,XLAMC,YLAML,RHCX,KHCY,ELAS) 

WKI  Tc.  IE,  3r  : ) NFC  P C T , KF  0 1 R , L SI  , L SN 
WRITE  (fc,420  RR,  XL  , XH  ,ElAS  ,XNI  ,UD  , EXXP 
WRITE (6*573) XcAMD,YLAMC,EXX,EYY, RHCX,RHOY 


73 


r 


CfllL  1PPERF 

Q +4*  + + * + ******  + 

WRITE  < t, 5^  8) 

DO  85  1<=1»K1 
LK=  IK  - 1 

WS  1 Tc  (6,:>lr)LK 
WRITE  ( £ , 5 2 r ) 

XX  = P . 

CO  85  1 1 = 1 » N EQP 0 T 

WRITE ( 6 , 5 G 9 1 I1,XX,WZ(I1, IK) » W'  2 F ( 1 1 » I K ) »WZPP(I1,IK) 

X X=  X X »lELT  a 
85  CCNTIME 

IF  (LfMM.KE.il  GO  TO  39 
WRITE (6. 5 CL ) DELT A, AL1 »GA1 , AL2,QT  2 ,GA2 
WRITF(£,5C!)Pll,Wl2*H22,G^l,Qi?,Q£2 
WRITE  { t , 5 ? 2 ) Ul 1 • Cl2»U22,GB2,DE3»  OB4 
WRITE(t»5r2)0Ll*DL2,DL2»DL4,Cl5 
WRITE  ( 6 , 5 3 4 ) OA  1 , CA2*DA3»DA4 
39  CONTINUE 

DO  63  11  = 1 * N E Q P 0 T 
DO  64  „i=l,Kl 

wm< ii , ji j = : . 

FTM(Il.Jl) 

WPP < II  , Jl > =3  . 
u 4 CONTINUE 

CO  65  j 1 = 1 . K 2 
FM ( II, jl ) =C  . 

XFMI 1 ,J1 ) =0  . 

FPF (11 ,J1 > =9. 

65  CONTINUE 
62  CONTINUE 

READ(5,2In)DUND, XNXX,XFRLE,DNXX, ACCUR, RII 
1RR  = R ] I 

IF  ( IRR  ,c.Q  . ) IRR=  1 
10LNu=  LUND 
CN  X = CN  XX 
x°«ES=XPR£E 
XFN  X=  X NXX 

C XFNX=AXIAl  COPFRESSICN,  XPRES=HYCR03  T AT  i C PRESSURE,  XNX  = .ITHE 
C uR  XPRES  ^CCCkUING  TC  IOLNC 
UO  TC  l 71  , 72 , 7.3  ) , ICLND 
71  WRITE (t, 511) xP RES, XF N X , DN X , A C CUR 
X N X = X N XX 
XN11=XNXX 
GG  TC  74 

«.J  1 1 (fc  , -j  1 2 ) XFNX  ,XFRES  ,ONX  , AC  CUR 
^f.A=XFFc.L 
«M  1=  * F RE  S 
44  TO  74 

- r (6,51-1 XNXX,XFRES,ONX,ACCUR 

i • • »N  « X 

' < • X F R t.  S 


■S 


XFNX 


# • is 


h^AU(5.1J^)NNN,LNNN.ILNW 
NWA VC=NNN 

**#♦»»♦»*♦*♦«*♦****» 

CALL  C CFF  NN ( NW A V c ) 

*»♦**«*♦***•*♦**♦**** 

HR  i Tr  < c , 5 " 5)  NNN 
IF ( LPRINT . NE  .1 ) GO  TO  49 
WRITE (6*50c)C1»C2*C3»C4,G5»C6 
WRITE (6, bC 7)C7,Cfl,C9,C10,CIl  ,C12 
49  CONTINUE 
1LR  = G 
L I C C N = 1 
IPO  T T = 0 

CALL  SECOND  ( TIFI1  ) 
nR I T c ( 6 , 7 9 3) TIPI 
T IM2=  T IM1 
TIF4=TIM2 
IINN=? 

555  L N - 1 

IF(IUlNO.EC.1)XFNX=XNX 

IF  (IClND.EO.E.  0 R^j  r-L  N C . EG . 3 ) XPRES-XNX 

IF(IULNJ.EC.:»XFNX=TlMDx*XFRES 

IDE  1 = ILOET 

CALL  F CTL  Fo ( IDET  ,NRHS» PAXN, AP.EP , CP,GP,PR,XP,CC»HT,Tl,Vi,MAX2, 
1 1 X F M * C FTP*  XFNX»LN'«NJ»NW*NF  ) 

IF ( LPk 1NT  .HE  .1 ) GO  TO  lni 
CALL  SFCONr(TTP3) 

T I P 1 = T Ir13 - T IP2 
TIM2=  TIM  3 
T I P 4 = T IM  ^ 

Wki.Tr.  < 6 » 2C1INHAVE*  XFNX  , XPKC  S » TI M 1 
::i  CALL  T RANSF ( WF , XWF ,FF, XFF.NW ,NF, 2 ,Tl , M A XN  , 1 , LP  fi  i NT ) 

-*•♦4  luET  = ILOET 

IF(IlLF0.E0.1)XFNX=XNX 

IF  UCLNU  . EC.  2.GR.IoLND  .EQ.  3)  XFRE  S = XNX 

IF  CILLNU . EC.  7>  XFNX=TLMDX*XPRES 

IP AX=  1 

WP AX="  . 

ITER-C 

co  ir  <.  ji=i,ki 

1C  2 WMA  X=*PA  X + WP  (1 , J 1) 

CO  lr2  I 1 = 2 ♦ NEOPCT 

wpp  = C . 

DO  104  J1=1,K1 
1*  4 h M M = W tv  P ♦ WM  (II,  J1  ) 

IF (AfcS (wMP)  .LE.AES(WMAX))  GO  TO  1C3 
W”AX=NrM 
IM A X=  1 1 
i: 2 CONTIN  Lr 
J « P A X = 1 

wWP.  (1  ) = WM  ( I PAX  , 1 ) 

AW»*P=WWM(1> 

IF(Kl.cQ.l)  GO  TC  1351 

DO  x r 5 J 1 = 2 t K1 

WWl‘ ( J 1 )=HP ( IPAX , J1 ) 

IF  < *9S  (WWM  ( Jl)  ) . LE  ,A3S  (AWWrt)  I GO  TO  105 
A W W P = W Wil  (Jl  ) 


J A i’  W X - v 1 
lr  5 CONTINUE 
1151  JF  N A X = 1 

FFM  (1 ) - FM  ( It' AX,  1 ) 

AFT  N=FFM( 1) 

IF  ( K2 . EQ  • 1 ) GO  TC  333 
CO  irt  J1=2,K2 
F F f ( J 1 ) = F N l I F A X , Jl  ) 

IF  (ABE  IFFMJU)  .LE.ABS(AFFM)  ) GJ  TO  ICE 
AFFF-FFM  (J1  ) 

JFM AX  = Jl 
106  CONTINUE 
333  LN=  2 

IT-_K  = ITEP  + 1 

IF  (II  Eh. LE .1.)  GC  TO  113 
*P 1TE  (6, 1 IN)  ITER 
GO  TC  999  9 

113  CALL  FCTERS (IOET  ,NRHS , MAXN , AF  ,bP , CP,GF,PR,XP,CC,MT,T1,V1.MAX2» 
1IXPM, CETM.XFNX ,LN,NJ  .NN,NF ) 

IF (LFR1NT.NE.1)  GO  TO  111 
CALL  SECONU(TIN3) 

T I M 1= 1 IM3-TIF2 
TIM2=TIM3 

WRITE  (6*  lid)  ITE'-’  , N N A V F.  , XF NX , X FRE 5 , T I Ml 
111  CALL  TFANSF  ( NF  , X HF  , F F , XFF ,NW,NF, 2 ,T  1 , M A XN  1 1 , LP  ft  I N T ) 

CO  116  J 1 = 1 » Kl 

IF (WM  ( IMAX  ,J1>  .NE.  r. ) GO  TO  57 
WCCN (1 TE  P , Jl ) = 0 . 

GO  TO  115 
57  CONTINUE 

WCON(ITER,Jl)=ABS(  <WM  ( IMAX  ,J1)  - WN’MtJl ) )/WM  ( IMAX,  Jl)  ) 

.1 5 CONTINUE 

WCH=WCCN<ITEK*JWMAX) 

IWH-JWNAX 

oo  lie  ji=:,k2 

IF (FM  (IMAX,  Jl)  .NE.  : . ) GO  TC  5e 
FCCN  < I TER, Jl>  = ° . 

GO  TC  116 
53  CONTINUE 

FCCN(I1EP,J1)=ABS((FM(IMAX,J1)-FFM(J1))/FM(IMAX,J1)) 

lie  cgminle 

FCH=FCCN ( ITER, JFMAX) 

IFH  = JF  FAX 

IF (LFP 1NT  .NE .1)  GO  TO  117 
WRIT-:  (6,  llo)  ITER  ,WCM,FCH 
W c I T f (6,119)  (Jl, WCCN(ITER,J1) tJl-l*Kl> 

WRIT'  ( c , 119)  ( J 1 t FCCN ( I TER , Jl ) , J 1 = 1 » K 2 ) 

117  IF  (WO  .GT  .ECCNV)  GC  TO  194 
If  CFCH.GF  .ECCNV)  GC  TO  194 
GO  TO  195 

194  IF  ( ITEfi.LF.2)  GO  TC  196 

IF (WCCN ( I Tt R * IWH  ) . GT. WCON( ITEk-l , IWH) ) GO  TO  197 
IF <F^UN(  lit R,1FH  ) . GT  .FCON ( ITER-1 , IFH) ) GO  TO  197 

GC  TO  196 

197  IF  (XNX.NE.XM1)  GO  TO  198 
WRITE (6. 991 ) XNX 
GC  TC  9999 

I ( 

[ I 
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196  DO  151  J 1 = 1 * K1 
131  WWM (J1 )=  WM ( INAX , J1 ) 

.DO  132  J 1 = 1 , K2 
152  FFN ( J 1 ) = F N (INAX, Jl) 

GO  TO  533 

195  IF(IOLNO.EQ.l)XFKX=XNX 

IF (1DLN0.E0.2.UR.ICLND.EQ.3) XFhES=XNX 

IFdO-.NO.  EQ.3)  XFNX=TLMCX*  XPRES 

CALL  FCTSN(FCT,STRY,STRA,..,1,1,XFNX) 

CALL  SEC0NDITIM3) 

T I Ml  = T 1M  3 - T I M4 
TIM2=T  IM3 
TIM4=TIM3 

WRIT-C6,241>XFNX,XFRES,NWAVE,iTER,TIMl 
WRITE (6*242)POT  , STRY,STRA 
IFdUcT.EC.l)  4 RITE (6, 243) DETM,I XFM 
IF (LMLl.NE.1)  GO  TC  476 

CAlL  TRANSF  (KF  , X VsF  ,FF  , XFF  , NW  , NF , 2 , T1 , N AXN  , 2 , 3) 

476  CONTINLE 

IF  (LNNN.  EC  .2  . AND  .LICON  .NE.  10  ) GO  TO  666 
IF (LICCN.NE.1DI  GO  TO  629 
ILR  = 1 
GO  TO  777 
629  XNX1-XNX 

1 1 NN= 1 1NN  + 1 

IF (IINN.LE.IRR)GC  TO  721 
WRITE ( 6,7  2 2 ) 1INN 
GO  TC  9995 
721  CONTINLE 

XNX=XNX+DNX 

IF (TXNX.GT  ,XNX)  GO  TO  244 
DN  X = L)N  X/2 • 

XNX=XNX-GNX 
ACN=CNX*1C( / XN  X 
IINN=IINN-1 

IF ( ADF ,GT . ACCUR ) GC  TO  244 

XNX=XNX1 

GO  TC  «19 

244  IF (INCN.EC.l)  GO  TC  555 
REWIND  2r 

WkITE  (29  ) ( ( W M ( 1 1 ,J1)  ,J1  = 1»K1)  ,Il  = l,NEQFOT)  , <C£T(1til  ,J1)  ,J1  = 1»  Ki) 
1,I1  = 1,NEQPCT) , ( ( HTF(I1,J^)  , J 1 = 1 , Kl) , 11  = 1, NEGRO  T ) , ( ( FM( 1 1 , Ji) 

2,  Jl  = l , K2  ) , 11=1 , NEQFOT-)  , ( (XFM  (II,  J 1 ) , Jl  = l,K2>  , 1 1=1  ,NEQPOT)  , 
2<<FMPdl,Ji),Jl  = l,K2>, 11=1  ,NECPQT  ) 

GO  TO  444 

19  ft  IF (LILCN.Nl.1C)  GO  TO  429 
I L R = L' 

GO  TO  777 

425  IF  (LNf  N.  rO.J)  GO  TO  24<= 

IPCTT  = IpO  T T + l 
IF (IPCTT.GT.l)  GC  TO  249 
POTT=PCT 
F X N X=  X N X - C N X 
249  ADN=L)NX*lCr  ,/XNX 

WR  1 Tc  (l, 545) N WAVE, xN X 
T X N X=  X N X 
XfJX  = XN  X-ONX 


1 1 NN=  C 

IF  (AON.LE.ACCUE)  GO  TO  819  , 

CN  X = 0 N X/2  . 

XNX=XNX+CNX 

IF ( INCN. EC . U GO  TO  935 
REWIND  2 0 

PE  AD  <2  C)  ( CWn  ( 1 1,  Ji)  , ul=l,Kl), 11=1, NEC POT)  , t (ETfl  CI1.J1J 
1,1 1 = 1 ,NEQFOT) , < < WFF (II , JI) , Ji=l, Ki> , 11  = 1, NEUPO  T)  , ( < FM ( 1 1,  Ji) 

2,Jl  = i,K2)  ,Il  = l,NEQFOT)  , ((XFM(Il,Jl>,Jl=l,K2),Il  = l,i»E0F0T)  , 

2 ( (FFF ( II, JI ) , Jl=l, K2)  , 11=1 ,NECPOT> 

GO  TO  444 

319  IF(IULND.EQ.l)  XFNX=XNX 

IF ( IGlND.EG.2.0R .IOLND.EO. 3) XFRES=XNX 
IF  (IGLND,FC.3)XFNX=TLMGX*XPRFS 
WRITE  (6,  785)NwAVE,  XFNX,  XPPES,  FCT  , STft Y, STRA 
L C A LCU  L A r ICN  OF  CRITICAL  WAVE  NUFEER 
56c  WR1TE<6*29) 

IF ( LNNN. EG  . 0 ) GO  TO  9959 
WRITE  (6,584) 

ILR  = 1 

NWPRj.N  = NWAVE 

NWAVL=  NWAVE  + l 

IF (LNNN.EQ.2)  FOTT=PCT 

FG  T F I N =PO  T T 

NM I N=  A NN 

INWAVE  =0 

I S TOP=  G 

19  = 0 

FC  T=PC  T T 

77  7 IF ( IDLNO. EC. 1) XFNX-XNX 

IF (IClND.EC • 2 » OR .ICLNr.EO.3) XFRE  S = XNX 
IF ( IDLNO. EG. J) XF  NX=TLMCX*XPRES 
IFTILF.EQ.I)  GO  TO  778 
WRITE  (6,582 )NWAVE, XFNX, XPRES 
GO  TO  5999 
776  INrtAVE=INWAVE+l 
1 1 N N = C 

IF (LNNN.EC.2JPXN >=XNX 

X N X = P X N X 

LICCN=in 

IF  ( INWAVE . EC .1)  GO  TO  291 

NWPRIN  =NW A VE 

IF (19. CE. 1 ) GO  TC  6 9 4 

IF  (PCIMN.Lc.PCT  ) CO  TC  139 

FOT  FIN=P0T 

NMIN-NFAVE 

N Wm  Vl  = N W A VE  ♦ 1 

GO  TO  291 

139  IF  (NWAVE. LE.NNN*  1)  GO  TO  549 
1ST  OF  = 1 
GO  TC  133 
r>4  9 15  = 19+1 

IF  ( 19  . GT . 1 ) GO  TO  694 
NWAVE  = MN-1 
GO  TO  291 

694  IF  (PL  tFIN.L5.FCT  ) GO  TC  095 
PGTMIN=POT 


1 


NMIN=NRAVE 
NWAVE=NWAVE-1 
GO  TO  291 

u95  ISJOP-1 

GO  TO  185 

39 1 CALL  CCEFNN(NWAVE) 

CALL  SECONC(!IM2> 

18  5 rtRITE  (c,581)NWFRIN,PXNX,POT » T IM2 
IF  (NWA  \E  . Lt . l ) GO  TO  9599 
IF (I3TCP.NE.1)  GC  TO  758 
WRITE <6,789)XNX,FOTMIN,NMIN 
GO  TO  5999 

79e  IF ( IN WAVE  • G T . I L N W ) GO  TO  9995 


GO  TO  555 

5999  kFAO ( 5 , lr ° > LCONTN 

IF  (LCCMN.EG.l)  GO  TO  1111 
2 2 FORMAT  ( 1 H 1 ) 

10  FORMAT  (1HC  ,5A8) 

bo  Format (//25h  beginning  of  next  case  ///) 

ICO  FORMAT  <19  16) 

20  9 FORMAT  (6E12.A) 

3CG  FORMA!  {// ,2X, "NO  . OF  P G IN T S=" , I 3 , 2X , "KFOU R=",I 8 , 2 X , "oO L NO . LON  OF 
1FGINT  l-",ie,2X,"BCUNn.CON  CF  POINT  NECPOT=“,ie> 
h0  0 FORMAT  (//,*_X,"R  = ",El2»A»2X»"XL=,*»El2.  A , 2 X , "XH= " , E 12  . a,  2 X, 

1"EL AS-", E 12.  A, 2 X ,”XNI  = ", El2. t,2X , "JD  = ” ,El2. A, 2x,"EXX^  = " ,El2 . A > 

50  8 FORMAT  (// ,2X, "THE  IMPFRFECTICN  FORM  IN  AXIAL  DIRECTION  I S’*/ » 

510  FORMA  I (// ,2X ,"THE  IMPERFECTION  FOR  CIRCUMFERENTIAL  WAVE  ",I6/) 

52  0 FORMA!  {//, AX,  “POINT" , 9X,"L  ENGTH" , lcX ,"WZ" , 1 AX , " W ZP" , 1 2X ,"WZPP"/ ) 

50  5 FORMAT  (I  1C  , AElt . 6) 

5C  u FORMA!  (//  , 2 X ,“C E IT A = " , E12. A , 2 X, " A L 1 =" , E 12 . A , 2X , "G Ai =" , E 12 . A , 2 X , 

1" A L 2= ",  El  2.  A ,?X, "0! 2 = ", El l. A,2X,  "GA2  = " , El  2. A) 

50  1 FORMAT  (//  , 2X,"Hi  1 = ",E12.L, 2X,"H. 2=", El  2 . A , 2x,"H22=" , - 12 . A, 2X, 

1"Q 1 1="  ,El 2 .A  ,2X , "0*2  = " ,Ei2  «A ,2X,  " C 2 2 = " ,El2 ,A) 

50  2 FORMAT  <//  , 2X  ,"C1 1 = ",E12. A, 2X ,"U12=", tl2. A , 2X,"C22=" ,-12 . A,  2X, 

1"0  b 2 = *' , F 1 2 .A  ,2X  , “OE3=",E12.  A ,2X,  "OeA=" , c 12 . A) 

5C  3 FORMAT  (//,  2 X,"OL  1=",  El  c . A , 2X  , "DL  2=" , E 1 2 . A , l X , " DL  3 = " , - i 2 . A , 2 X , 
1"ULA-",E12.A»2X,"DL5=",E12.A) 

5P  A FORMAT  (//  , l X ,"  C A 1=  " , E 1 2 . A , 2X  , "DA  2 ="  , E 1 2 . A , 2X  , "0  A 3 ="  , - 1 2 . A , 2 X , 

1"DAA  = " ,E12 . A ) 

,11  FORMA  1 (//,  2X ,"F OR  FIXEC  PRESSURE  OF  " , E 1 2 . A , 2 X , " T HE  INITIAL  AXIAl 
1LOAIJ  IS  " ,l12.  A , 2X,"THE  INCREMENT  OF  AXIAL  LOAD  IS  ",E12.A/ 
1lX,"ANC  THE  ACCURACY  (FERCENT)  OF  THE  AXIAL  LOAD  IS  ”,E12.A) 

512  FOEMAI  (//,2X,"FOK  FIXED  AXIAL  LOAD  OF  " , E 1 2 . A , 2 X , " ! HE  INITIAL  H YGR 
1CSTAIIC  PR  r SSUR  E I S" , E 1 2 . A/ 2 X , "T H t INCREMENT  OF  THE  HYCROSTATIC  PR 
2 E S S U F:  E I S"  , E 12  . A , 2 X , " AND  1 HE  ACCURAoY  (PcRCENT)  OF  THE  HYDROSTATIC 
? PRESSURE  1S",E12.A) 

50  5 FORMAT (// ,2X,"TH=  CIRCUMFERENTIAL  WAVE  NUMuER= ” , I 6 ) 

_»C E FORMAT  (/  / ,2X,“Cl=",r 12 ,A, 2X,”C2=",ri2. A,rX,"C3=”,E^  2. A, 2X, 

1 "C  A -“ » E 1 2 . A,2X,"C5=",£12. A,2X,"Co  = ",212.h> 

5o  7 F C-MA  T (//  , 2X ,"C7=" ,E12  .A, 2X,"C8= " ,E1l . A , 2X,"C9=", E12. A , 2X, 

1" C 1 l = " , E 1 2 . s , 2 X , " C 1 1 = " ,E12.A,2X,"C12=",-12.A) 

79  3 FORMAT  (//,  2X  ,"rL flPSED  T I M E=" , E12  . A , 2 X , "SE CONOS’V  ) 

cC  1 FORMA  I (// ,2X  ,"  IN  IT  iAL  SOLUTION  (FFlM  LINEAR).  FuK  N =" , I 8 , c X , " N X-" 
i,E12.  A,2X  ,"F=",E  1L  . A/2X,"TI.1E  COPUATICN  =",Elc.A,t  X , "S ECONO S" / 

22  X ,"  - - r=  = = r = = = s = sr  = ==  = = = = =.=  =*  = = = = = = = = = = = = ===  =*  = =.=«.=  = = ==  = == 

11  A FORMAT  (// ,2X, "END  CF  THIS  CASE  BECAUSE  ITER  GklAIEk  THAN  ", 


L 


18) 


112  FORMAT  <//,  2X  ,"SOLUl  ICN  FOR  ITERATION  " , 18 , 2 X , " N=" , i 8 , 2 X , "NX=" , 8 

1E12.4,ZX,"F=”»E12«4/2X,"TINE  COM  PUT AT ICN  =",E12.  4, 2 X, "SECONDS"/ 

2 = = = = ^ = ^ = = = = = = = = = = = = = = = = ==  = = = = = = = = = = = = = = = = •*) 

11 C FCrMAT  (// , 2X ,"ITER=", I 8,2X,"WCH=" ,El5.5 ,2X,"FCH=",-  lv. 5/1 
991  FORMAT (//,2X, "THE  SOLUTION  IS  NOT  CONVEKGEL)  AFAKENTLY  THE  INITIAL 
lLCuU  SHOULD  fcE  LESS  THAN  ",E12.4> 

241  FORMAT (// ,2X, "NONLINEAR  SOLUTION  FOR  N X=” , E 1 2 . 4 , 2 X , "P  = " , E 12 . 4 , 2 X , 

1"N  = " » I 8 , 2 X , " IS  OBTAINED  BY  " , Io , 2 X ,"  I T ERA  T IONS  *V  2X , ”T  I M E CcMPLJAII 
2CN=",E 12. 4 ,2X,"SECCNOS"/2X ."CKCKOKOKUKCKCKOKOKCKOkCKJKCKOKOKOKCKJK 
2lKCKCKCKU(<CKCKCKCKC<OKCKOKCKCwOKOKOKOKCKOKCk>OKOKUi\CKOKOKQKO<0“) 

242  FORMAT  <//, 2X, "POTENTIAL  ENERG Y=" , E 15 . 5 , 2 X , "END  SHORTENING  FOR  Y = ~. 

1 IS", cl5. 5, 2X, "AVERAGE  END  SPCRT E N ING=" , E 1 5 . 5 ) 

24  7 FORMAT  <// ,2X,"DETE KIN ANT=",E12.4 , 2X,”IXPM=", 1 10 ) 

76  5 FORMAT (//,2X, "CRITICAL  LOAD  FOR  N=" , I 8 , 2X , " I S N X = " , E . 2 . 4 , 2 X , 

1"P=",l12.4/2X,”PCTENTI AL  ENERG Y= " , El 5 . 5 , 2 X , "END  SHORTENING  FOR  Y=’ 

2.  =", t 15 . 5 ,2X," AVERAGE  END  SHORT EN IN G=" , c 1 5 . 5/ 2 X , "G KOKC KOKOKO KO KOK 
ICKCKCKCKOKCKCKOKCKCKOKCKOKOKCKCKOKOKUKCKOKOKOKCKCKCKjKCKOKUKOK") 

11 9 FORMAT  (5(18, cl  5. 5)  ) 

545  FORMAT (///2X,"FOk  N = " , 1 8 , 2 X , " THE  LOAD  ", El2.4,2X,"IS  OVER  THE  LIM 
1 I T PCTNT") 

57  2 FORMAT (// ,2X ,"XL AMC  = " ,E12. 4,2X," Y L AMO= " , E i 2 . 4 , 2X,"-_ X = " , £12 . 4 , 2X , 
1"EY=",c12.4,2X,"RHCX=”,E1l.4,2X,"RH0Y=",E12.4) 

•38  2 FORMAT  (//,  2X  , "FOR  N = " , 18 , 2 X , "NX=  " , Ei  2 . 4 , 2 X , "P=“  , El2 . 4 , 2 X ,"  T HE  SOLIT 
1IGN  IS  NOT  CCNVERGFO"/2X,"RRCBA9L Y THE  CRITICAL  LOAD  IS  SMALLER  TH 
2 AN  TH  r SE  LOADS") 

7o  9 FORMAT  (// ,2X, "FOR  =" , E 1 2 . 4 , 2X , " THE  HINIMUN  POTENTIAL  ENERGY  IS  " 

1 » El  5 . 5 ,2X  » "AND  THE  CRITICAL  RAVE  NUMBER  IS  ",I8> 

38  4 FORMAT  (/ / * 2X  , "D  E Tr  Rtt  1 N A T I ON  CF  THE  CRITICAL  CIRCUMFERENTIAL  WAVE 
1 N UNDER 

7,  *»»♦,»»»♦******••) 

512  FOkMAT (//,2X,"THE  RELATION  0 E T WE  EN  AXiAL  LOAD  AND  THE  FRESSURE  IS 
1 NX=",l12. 4, 2X, "MULTIPLY  LY  THE  P RESSU R E"  , c 12  . •♦/ 2 X , 

2"THE  IRCRlMENT  OF  THl  FRtSCURE  I S“ , c 1 c . 4/ 2 X , 

3"AND  THE  aCCCRACY  (PEF  CENT  ) CF  THE  PRESSURE  IS",E1Z.-,) 

5h 1 FORMAT  <//,2X,"WAVE  ND  =" , 1 1 ? , 2X , " LC AD  = ” , E 1 5 . 5 , 2 X , 

1"FCTEN  TIAL  £ NE  R G Y = " , l 1 E . 5 , 2X  ,"EL  AFSED  T IM  E = " * E 1 l . 4 * 2 X , " SECOND  S" ) 

722  FORMAT (//, 2X, "END  OF  TFIS  CASE  3ECAUSE  NUMBER  OF  lCAU  POINTS  IS 
IlFEAIER  THAN", 18) 

ENO 


i 


SUCiKOlllNF  I NPE  R F 

CCMMCN/FOURIR/KF  C U R » < F »K4,K3»K2»  K1 
GOMNCN/CIN1G/NEQFOT,KI  (5rC  ) 

COMMON /FACT  3/01 5 , XL,XH 

COM1CN/FIlFR/DELTA,AL1  , G A 1 , Ate  , :3  T 2 , G A2 
COMMCN/PRES2/WZ  ( 2 3*  ,5)  ,WZP(2!T.  »5>  , WZ FP  ( 2e  " , 5) 
DV  = 1. 

PI=4. ’AT  AN  (CV) 

AC- 1 . L *XH 
Al  = PI’ 7./4  . 

X X = C • 

00  10  I1=1,NEQF0T 
N Z ( 1 i * 2 ) =-AC*SIN  (Al’XX) 
WZP(11,2»=-AC’A1*CCS(A1*XX) 
WZPP<Ii,2)=AC’Al’Al*SIN<Al*XX) 

W Z ( 1 1 , 1 ) = " . 

KZP(I1 » 1 ) = c . 

VsZFF(Il»  1)  = r. 

IF (Ki  .LE . 2 )GC  TO  50 
00  11  J1=3,K1 
WZdl.JJ  = C . 

wzp(ii  ,jd  = ?. 

WZFF(I1,J1)=C. 

11  CONTINUE 
5 0 CONTI  ME 

XX  = XX  +CEL  T A 
10  CONTINUE 
RETURN 
END 


SLMRULlINE  TrANSF  ( NF , XKF,  FF  , XFF  , Nfc.NF  , NRF  , T 1 , M A XN  , 1 Dc R , 1PRK)  . ~ 

CCMMCN/PRES1/NM  ( 2 ? , 5)  , E T M ( 2 C , 5 ) , WMP(2:3,5)  A C 

CCMM0N/PRES3/FM  8J  ,XFK  ( 2 G T , 0 ) , FM  P ( c ' r.  , 8 ) 

COPrlCN/FILFS/b’-.L  TA  ,AL1,GA1  » A L 2 » B T 2 » G A 2 
CONMCN/CCIS</I2i ( 5 ” 1 ) , 122  (501) 

common /Ci  mg/nfqfot.  mi  <5:  r ) 

CGHMCN/FOURIR/KF  CUP  » KB  , KA , K3  , K2 , K1 

DIMENSION  N F (NKF  , N W ) « X K F (NKF  » N W ) , FF  (NRF.NF  ) , XFF  (NRF  ,NF  ) (MAXN) 

IF  (IPrtR.cG.3i  GO  TG  272 
DO  1C  1 1 = 1 , N EQP  0 T 
NL  = M 1 (II  ) 

CALL  Ft  ADM  S(  22.T1.NL,  ID 
IF  (Il.NE.l)  GO  TC  175 
CO  11  J1=1,K1 
WF(i.Ji) = T 1 (Jl) 

WMCll.wl ) =T1 (J1+K6) 

XWF (1 , Jl)  = II (Jl  + K5-1 ) 

ETM(Ii.Jl) =T1(J1+K3+K6-1) 

11  CONTINLE 
00  12  J 1 = 1 , K 2 
FF(l.Jl) = T 1 (Jl  + K 1) 

FP(ll.Jl) = T 1 (J1  + K1+K6) 

XFF (1 , .11 =T1 (J1+  KA ) 

XFM ( II , Jl) = T 1( Jl +KA+K6 ) 

12  CONTINLE  - 
GO  TC  in 

175  CO  12  wl=l,Kl 

WM(I1,^1)=T1(J1> 

ETP (II ,J1) = T i ( Jl +K3-1) 

1 J CONTINLE 

DO  1 A oi=l.K2 
F P ( 1 1 , Jl  ) = T1  (Jl+Kl) 

XFM  (1 1 ,J1)  = T1  (Jl  *KhI 
1A  CONTINLE 

IF (Il.NE.NLGrOT ) GC  TO  1C 
DO  15  J 1 = 1 , K1 
WF (2  , J 1)  = T1  ( Jl +K6) 

XWf(2,jl)=Tl(Jl*Ki*K6-l) 

15  CONTINLE 

DO  lc  J 1 = 1 , K 2 
FF (2.J1) = T 1 (J1+K1+K6) 

XF  F ( 2 , Jl  ) = T 1 (J  1 ♦ Kh+K6) 

IE  CONTINLE 
IP  CONTINLE 

IF  ( IU_K.NE.1I  GO  TC  275 
NEt.P-NrOPCT-1 
DO  11  = 2, NEQP 
DO  21  „1=1,K1 

NMF(I1,J1)=AL1*WMI1-1,J1)  ♦GA1*WM  (Il  + l.Jl) 

21  CONTINLE 
DO  22  o 1 = 1 , K 2 

FMP(Il,Jl)=ALl*FP(il-l,Jl)4GAl*FM  (II +1  , Jl) 

22  CONTINLE 
2"  CONTINLE 

DO  2o  J 1 = 1 , K 1 

.lMF(l,<,i)=ALi"NF(l,Jl)4GAi*WM(2,Ji) 
KPF(NECPGT,Jl)rALl*W/i(NEQP»J3>+GA.i*WF(2,Jl) 


23  CONTINUE 

CO  24  v 1 = 1 , K 2 

FMF(i»Jl)  = AL 1*  F F (1  iJ.  ) +GA1*FN(2*J1) 

FMP  IN:  CPOT  , J 1)  = All *FM (NEQP, Jxl +3Ai*FF<2,Jl> 

24  CONTIMF 

275  IF  (IFRk.NE,l)RETlRN 
278  CONTINUE 
J1  = C 

WRITE  IE, 4 3 . ) J1 
WRITE  (F,  50  0 
X X = 0 . 

WRITE (6, bO  > WK (1 ,1) , XWF (1 , 1) 

DO  48  11=1*  NrQP 0 T 

WRITE  (fc ,509)  11, XX, WM ( 11,1) , WKP<I 1 ,1 > ,ETM< 11,1) 

XX=XX+ CEUTA 

48  COuTIME 

WRITE  <6,60 C) WF (2,1)*XWF(2*1) 

DO  49  ol=l,KFCCR 
WRITE  < 6 » 4 I C ) Jl 
WRITE  <6, 52C  > 

WRITE  (6, 7 P C ) WF (1  * J 1* 1 ) ,XWF ( 1 , J 1* l),FFCx,Jl) ,XFF (1 ,J1) 

CO  51  11= 1 , N EOFO  T 

WRITE  (£,or9)Il»WI'(Il,Jl*l)  ,WNMIi»Jl*i)  ,£TM(II,Jl+i),FM  (II, Jl)  * 
1FMF  (II  ,J1  ) , X F M ( 1 1 , Jl) 

51  CGNTINLE 

WRITE (6, 7 0r ) WF (2  * J A -H  ) ,XWF  (2 ,Ji  + 1 ) ,FF  ( 2 , J 1 ) ,XFF (2 , J 1) 

49  CONTI  ME 

DO  52  wl=Kl,K2 
WRITE  < C » 4 C O Jl 
WRITE  ( 6 , 5 1 C ) 

WRlTu(e,oCi.)FF<l,Jl),XFF(l,Jl) 

DC  53  11  = 1 ,NEQP0T 

WRITE  <c,  7 '9)  I1,FM  11,  Jl)  , FC.P  (11,  Ji)  , XFM<I1 , Jl) 

5 7 CONTINUE 

WRITE  ( t , 8 r •.  )FF(2  ,J1)  ,XFF(2,J1) 

52  CONTINUE 

40  0 FORMAT  (// , 2 X ,"  I N TE  RM  ID  I AT  RESULTS  FOR  KFOUR=”,  18/2X  ,••»********,■»♦♦ 
5J  V,  FORMAT  (//,  2Xf"FOlNT",  4X,"l'ENGTH-  , i''Xf##W»  l4X,"KP",.3Xt"fcPP"t 

2**»4.*.f#*4#.**4HMA*¥»**MMMMHMH***#****MMMMMM*'*4MMMF*¥*4MM^**+************ 

5u 9 FORMAT  (18,  lie. 4,  2£15. 6) 

6C E FORMAT  (//2nh  FICTIVE  POINT  El 5 . 6 , 1 5 X , E 15 . 6// ) 

oC 9 FORMAT  <18 , 12X,6E 15 .6) 

7(i,  FORMAT  (//2CH  FICTIVE  FJINT  E15  . t , 15  X , 2ci  5 . 6 , - 5X  , E 15  . o/ / ) 

80  ^ FORMAT  <//65P  FICTIVE  POINT 

1 E15.6»15X,E13.6/7 ) 

7L  9 FOF.MmT  (Id,57X,3E15.6) 

RETURN 

END 


SUE  ROUTINE  A DUG  < 1EC,  Ml  ,CF,  L>F  ,AF,  bF  ,NPHS,  XNXX,  LN  ,NJ,  NW  ,NF> 
COMMGN/CINTG/NEQFOT,MI  (50C ) 4 £- 

COMMCNXJCUNC/LS1  ,LSN 

COMMCN/FICFR/DEl TA ,ALi ,GAi , AL2,9  T2.GA2 
UCMMCN/FOLPIR/KFCUF, KE ,«L,K3 ,K2, K1 

CIMENSION  AF  (Ml  , Ml ) ,i3F  (Ml, Ml)  ,UF  (Mi, Mi), GF< Ml, NRHS) 

C NEGPUT -NFOI  NT  (E  X C lUC  I NG  FIC1IVES  POINTS) 

C LSI  - K i N C OF  BOUNDARY  CONDITIONS  OF  POINT  1 
C LSN  -MMj  CF  EOUNOARY  CCNuiTICNS  OF  POINT  NP 
C NW  - MAXIMUM  K+l  FOR  DIMENSION  KM(NJ,NW) 

C NF  - MAXIMUM  2*K  FOR  CIMENSICN  FM ( NJ  ,NF ) 

IF (I5C.GT.1)  GO  TO  10 

CALL  *STG(FF,CF,AF,GF,1,XNXX,M1, N J , N W , NF , L N , NPHS ) 

DO  2 11=1, K6 

GF  < I1+K6 , NRHS) =GF < II , NRHS ) 

DO  2 J 1 = 1 , K 6 

9F(IltK6,Ji>=AL2«EF(Il,Jl)+ALl*GF(Il,Jl) 

BF<Il  + Kb,Jl  + K6)  = ET2*BF<Il,JiMAF(Il,Jl) 

OF  (II,  Jl*KV*'=bA2*BF(Il  ,Jj.)+GAi*CF(I1,J1) 

2 CONTINUE 

CAlL  ECUNCk ( AF,CF,GF, 1 ,XMXX, LSI, M 1 , N J , NW , NF , LN  , NRHS > 

CO  3 11  = 1,  K6 
DO  3 J 1 = 1 , K6 
OF  < II  ,wl) = A L 1*  A F (II, Jl) 

AF (I1  + K6,J1)=EF(  1 1 , J 1 + K 6 ) 
cF  (I  1 ,jKK6)  =CF  (II  ,J1) 

AF ( II , Jl) =G Al*  AF  (II, Jl  ) 

3 CONTINLE 
RE 1 URN 

IF  IF ( IEC  ,GT  .2)  GO  TO  2C 

CALL  R STG (AF, CF, cF, GF, 2, XNXX, Ml, NJ,UW.NF,LN, NRHS) 

DO  4 11=1, K6 
CO  4 Ji=l,K6 

DF(Il,Jl)=DF(Il,Jl)+dT2*AF(Ii,Jl) 

CF  (II  ,J1+K6)=AL2*AF  ( j.  1 , Jl ) * A L 1*  D F (II, Jl) 
AF(I1,J1)=GA2»AF(I1,J1)+GA1*CF(I1,J1) 

CF ( II , Jl ) =r  . 

A CONTINLE 
RET  URN 

22  IF ( IEU  .GE.NEGPUT-1)  GO  TO  2r 
J ° = I F.  C 

CALL  F S TG  (AF,CF , EF, GF ,JP, XNXX , Ml  , N J , NW , NF , LN , N RHS ) 

CO  b 11=1,  K6 
DO  5 J 1=  1 , Kc 

E F ( 1 1 , vi  1 ) = E F (II  , Jl  ) ♦ B T 2 * A F (II, Jl  ) 

TEMP  = u A2*AF  ( II , J 1) 4GA1*CF ( 11 , Jl) 

CF(Il,*.l)  = AL2*AF(Il,JlM-ALi*CF(Il,Jl> 

AF ( II , Jl ) =TEMP 
5 CONTINLE 
RETURN 

3 C 1 F ( I E L . Eu  . NECPO  T ) GO  TC  4C 
JP=  I EG 

CALL  SSTG(AF,CF,eF,GF,JP, XNXX, Ml, NJ,NW,NF,LN, NRHS) 

00  6 11  = 1 ,K6 

DO  6 j 1=  1 , K6 

CF  <11  ,J1)  = b F (L  , Jl)  ♦ 6 T 2*  A F ( I , Jl  ) 
rEM-  = LA2’AF  (II, Jl)  ♦ 3 A 1 * C F (II, ol) 


CF  (11 , 01  > =AL2*  AF  <11,  J1  ) ♦ALlMFtl  1 ,J1>  4) 

AF (II, J 1 ) = TEMP 
AF ( II » -1 f K6>  =C . 
t CON  T 1 ME 
RETURN 
40  JP  = IEC 

CALL  RSTO  ( AF  ,CF,  cF  ,0F,  JP,  XNXX,IU  , NJ,NW  , MF  ,LN,NFHS) 

CO  7 11=1, K€ 

CF  (Il+K6,NRHS)=GF(I1,NRHS> 

00  7 J 1 = 1 , K 6 

GF (II * 01 ) = BF (II, J1 ) + BT  2*  AF  (II, J1 ) 

EF(ll,ulfK6)=GA2*AF(ll,Jl)+GAl*CF(Il,Jl) 

CF<I1*K6,J1)=AL2*AF<I1,J1)+AL1*CF(I1,J1> 

7 CONTINLE 

CALL  eCUNOR (OF  , AF ,GF, JF,XNXX,LSN , V.l ,NJ, NW,NF ,LN ,NRhS) 

CO  8 11=1, K6 
TEMF  = CF (II  ,NKHS) 

GF (11 » NRHS  ) = GF (I1*K6,NRHS) 

GF (I1+K6,NRHS) =TEMF 
00  8 «1=1,K6 

OF  <I1*K6,  J1  + K6)  =GA1*CF  (Ii,Jl) 

CF  (Ii*K£>,Jl)=ALl*CF(Il,Jl) 

CF  ( II  , J1 ) =EF  (11+  K(j  , Jl) 

GF (11  + K6, J1)=AF ( II, Jl) 

8 CO  NT  I ME 
RETURN 
ENO 


FUNCTION  ALS3  (I,  J ,L  .iliJF.Nc.  »N3,N4  , LLI 
C N4=l  FCR  C l JP , I ♦ J ) OR  BUP.I-JI 

C N 4 = 2 FCR  G(JF,J) 

C L L = 1 FOR  W LL  = 2 FOR  F 

COMHCN/FOURlN/KFCUF,Kb,K4,K3,K2,Kl 
DIMENSION  b(N2,N3) 

IF<L.ol.3>  gc  TO  1 c 

11  = 14-  J 

12  = 11 

GO  TO  20 
10  11=IAES(I-J) 

12  = 11 

20  IF(N4.EQ.i>  GO  TC  120 
12  = J 

GO  TO  100 

12  j IF  ( 1 1 . IE  • K FOUR)  GO  TO  lOf* 

A L B = r . 

RETURN 

10C  IF (L.LF.3)  GO  TO  110 
t T A = 1 . 

IF < I. EC. J I E T A = 0 . 

110  GO  TC  (11,  12,13,  14,1-3,16)  ,L 

11  K1=I1**2 
GO  TO  17 

12  R1=J**2 
GO  TC  17 

12  fil=2.*Il*J 
GO  TO  17 

14  Ri=(<i.-ETA>*I1**2 
GO  TO  17 

.5  Rl= (2 . -c  T A ) *J**Z 
GO  TO  17 

lfc  IF ( I-w .LT  . C ) E T A =- 1 . 

Ri=-2.*ETA*J*I1 
17  IF  (LL  .EQ.  1 ) 12=12 -»1 
ALE=Pi*H (JF,I2> 

RE  TUrN 
END 


SUiiKOU  TINE  RoTGtfi, b,T,G,  JF,  XNXX,  Ml , N J, NW , NF , LN  ,NRHS ) j, 

C LN=1  FCR  LINEAR  LN  = 2 FOR  NONLINEAR 

C NJ  NAXIilLM  PCIMS  IN  MERIDIONAL  DIRECTION 

L.  Nh  MAX  IN  UN  KFOURM  l,F-  li  A X IMU  M 2.*KFOUR 

C X N X X mXIAL  COMPRESSION  LOAD  J3-  THE  POINT'  IN  MERIDIONAL  DIRECTION 
C 0 M .1 C N / F OLRIR/KF  CL'E  , KE  , KA  , K 5 , K2  , K 1 

COMMON /GLCM/RR,DD,Hll,H12,H22,ulltC12,C22, 011,012,022 
CCNMCN/FACTCR/Cl,C2,C3,C4,ClJ,C6,C7,C8,C9,C13,Cil,Ci2 
CCMiTON/PPcSl/WM  ( 29n,5>  , ETM  (2Gi,5  ) ,WMP(20r  ,5) 

C0MM0N/PRFS2/WZ  ( 2 9 r , 5 ) ,WZP  <20l ,5 ) ,WZPP(2rC,5) 

C0MMGN/PRES3/FN ( 22  ° , o)  ,XFN  <2GL , 8 ) , F M P < 2 ^ U ,8) 

COMMCN/XXLCAC/XPRES 

01  MENTION  F(MtMl)  ,SCM1,H1)  ,T  (Ml  , M x ) » G ( M 1 , NRHS  ) 

C K£-6*KF0UR+2 

C XA=4*XFCLR+2 
C K3=3*KFCLR+2 
C K 2=  2 * KF  C L R 

C K1=KF0LRU 
C C9=(NNN/RR) 4 4 2 
C C1=DD4P11 

C C2=2 .* CC4H124C9 

C C3=2.4C124C9 

C C4=UD4P224C9442 

C C5=d11/D11*C9 

C C6=:./  (RR40ll>4C9 

C C7=C9442/ (2. 4C11> 

C Co=G..2vC94  42 

L C1E-C9/2. 

C Cli  = ?.''L9*ni2 

C Cl2~-Ccc4C9442 

DO  1 Il=i,Kc 
G < 11, NRHS) =0 . 

UO  1 Jl=l  ,K6 
R(I1,J1)=C. 

S ( 1 1 * w 1 ) = 0 . 

T < T 1 » *.  1 ) =2  . 

1 CONTINLE 
01  = 9 

DO  2 Ii=K3,K6 
T < 1 1 , 1 1 ) = 1 . 

Jl- J1 ♦ 1 
R(I1,J1)=-1. 

2 CONTINLE 

C f-CCIL  IOR  ILfl  EOLATION  FOR  I = C 
R(l,K3>=CL4Hll+Qli442/Dll 
T (l«Ki)=2»*Dll/ (RR*C11)+XNXX 
T (1»1)=1./(RR4424Q11) 

G ( 1,  NR  PS ) = - X NX  X 4 KZFP<JF,i  ) 

G < 1 » N N P S ) = G (1»NP. PS)  + XPRES 
CO  6 v'=l,KFCLR 
IS=J4 4 2 
M = J ♦ 1 

T<1,n-+J)=T  (1 » K3  ♦ J ) -Cr>/2  .*IS*NZ(JP,M) 

S ( 1,  J*  *>  = S (1  ,J  + 1 ) -C5* IS’WiP ( JP,M ) 

T C 1 ? w ♦ 1 ) =T  ( 1 * J ♦ 1 ) - IS/ 2 • 4 (LS’KZrP  ( J F • M ) + C64 W / < J F , M ) ) 

T(1»K4*J)=T (*,KAOI+Cl  4IS*WZ (JP , M> 
T<l,N.4j)=T(l,Ki.OM„l'*IS*WZPP<JP,M) 
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S(l,Ki«J)=SCl,KlfJH-2.*Cli*IS*WZP(JP,M) 

IF(LN.EQ.i)  GO  TC  6 

T(i,*3+J)  = T(l,i<3*J)-C5/2.¥lS*Wf'<JP,M)+Cl'>»IS*Ft'<JP,J) 
r ( 1 1)  = T (1  ,J  + 1 > - IS/ 2.*  <L5*LTM  ( JP,MMC6¥  W'KJP.M)  I 

1 ♦ C i ' * IS  * X F t*  ( JF  , J ) 

S ( 1,  J + l)  =S  (1 , J*1 1 -C5*IS*WNP<JP,,1)  +2.*Clr*I3*FMP  (JP,  J) 
T(1,K44J)=T(1,K4+JUC1C*IS¥WP(J»,P> 

I < 1 .K14JI  =T  (i,Kl+J)  D*Ib*ETf'  ( JP  ,M) 

5 < 1 ,K1  *J)  =S  <1»K1  fJ>  +2.*C1C  ♦IS*W>1P  <JP,M> 

G(l,t,hFS>  = G<I,NRHS)-C5/2.¥IS*  IW1  ( JP,M)  *ETM  (JP,KI  + Wl-iP  ( J F , M ) 4 * 
1-C6/4  . *IS*  (WMJP  »M  )**  2 > +C1C  ♦ ISM  wMlJP.M  *XFM  <JF,  J)  +ETM  ( JP»M> 
2FM  ( JF,J>  «-2.*hMP  ( JP,M)  *FMP  (uP,J)  I 
6 CONTIME 

ECUIHBR1LM  EOLATIONS  FOK  1=  1 , 2 ,,,,,,, KFCUR 

00  3 11=2, <1 
1=11-1 

I S = I * ¥ 2 
1S2=1S,¥2 

Till, 1)=-C6*1S*WZ(JP, III 
T{il,r>2)  = -G5¥IS¥RZ(JD,Il) 

T C II , I 1) = C4»  IS2 
Till,  Il+Kl-U  = -Ce¥IS2 
T(I1,I1*K2-1>=-C2*IS*XNXX 
T (II, I1*K4“1)=C3*IS”1./RP 
R( II, Ii*K3-l)=Cl 
Fill, I1+K4-1)=-Q11 
C 1 = C 7 * IS  ¥ k Z (JP, I 1) 

G < 1 1 , N FHS  ) =-XNXX¥W7FP  (Ja, II  ) 

IFILN.EQ.il  GO  TC  63 
T III,  1)  = T (11,1)  - C6*IS*V'M(JF,  II) 

1 (Il,Ko)=T  (Il,K3)-Ct7¥IS*WK(JF,Il) 

Till, Il)=T(Il, III- ISM  C5*cTM(JP,il+C6*RM(JP, III 

6 1 = C 7 ♦ IS*  ( WN  (JP , II ) +WZ  (Jp,  II ) > 

62=07/2 . *IS 

c.3  = EZM'1(JF,Il> 

G ( Ii, N F H S ) =G III, N R HS ) - C 5 * 1 3 * £ TH  ( J F , i ) *RM(JP,I1 ) 

1 -Cb*IS*WM < JF,  1)  MPi  JP,I1) 

6 C CO  4 J = 1 , KFCLR 
JS=J*M 
M = J M 

1 (Il»o+l)=T  (Il,J  + l)  + El¥JS¥KZ(v.F,M) 

IFiLN.EQ.il  GO  TC  4 
T(Il,j+l)=T(Il,J*l)4-El¥JS¥WmjF,i) 

T < 1 1 » 1 1 1 = T III,  II  )+£2*JS*HH(JFtM)  * i Wil  ( JP  , ’4)  *2 • * WZ  ( JP  » i'll  ) 

Gi  II , NFHS  1 =G  <11 , NRhS)  +E1/2  .*  JS*WM  (JP,M)  *♦  2 + E3*  jS¥Wf!  ( JP,M) 

1* (WM(JF,M)  + 2.*V*Z  (JF ,M>  ) 

4 COGTIME 
DO  5 J=1 , K2 

T(T1,Km  + J)=T(I1,K4+J)+C1,:*<ALB(I,J,1,WZ,jP,NJ,NW,1,1)* 
lALEiI,J,4,h2,JF,NJ,NW,l,l)l 

T ( ll.Kl+JI  = T ( I 1 * K 1 ♦ J ) + Cl"MALE<I *J*2,WZPP»JP»NJ»NVi,  1, l)  + 

1 ALtiil  , J * 5 ,WZPF,  4 P , N J , N W , 1 ,1)  ) 

3(Ii,Kl*Jl =S (II, K 1 ♦ J ) +C1T  * ( A L9 ( I ,J,3,WZP,JP,NJ,NW,1,1)  ♦ 

1 AL*JCI,J,fc,WZP,JF,NJ,NN,l,ll  ) 

IF  (LN.EQ  . 1 1 GO  TO  5 

T<ri,K4+J»  = T<Il,(<4*J>fCl"*<AL!3{I,J,l,Kr,JF,NJ,NW,l,l)  + 

1 ALOi  I ,J  , 4 ,WP,„P  ,NJ,NH  ,1,  1)  ) 
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T < 1 1 ,K1  + J ) = T (II , K1 *J) ♦ C1C*  (AU3  (I , J,2,lTM, JP,NJ ,NW,_  ,1)  ♦ 

1 A L 6 ( I ,J,5,ETM,JF,NJ,NU,l,i) » 
S<ll,KDJ)=S<Il,Kl*J)+ClD<ALd<I,J,i,WPP.Jp*NJ,NW,i,l)  + 

1 A L d ( I ,J,6,h(-;F,  JF,NJ,\'k,l,  i>  ) 

G ( II,  NPHS) =G (1 1 , NRPS)  +C1C  * ( ( ALd( I , J . 1 , KM. JP,NJ ,NW,1 , 1) ♦ 

1 ALO ( I ,J , 4 ,KK, JP ,NJ, NW,1, 1 ) ) *XFM ( JP, J)  M ALU (I , J ,2 ,ETM, JP,UJ,NK, 
21,  1 ) ♦ ALB<  I , J ,5,  E TM  , JP,  NJ,  NW,  1 ,1)  > ¥FM  < jF  , J >♦  (ALU  (I  , J , 3,  fc.MF,  JP,N  J, 

3 NW ,1, 1) + A LG (I* J ,6  ,WMP,JP,NJ, NW, 1 ,i> ) * FMP ( JP, Jl) 

IJ1=I«J 

IF ( IJ1  .GT  . KFCUR)  GC  TO  8r 

I (II  , I w I + I » = T (11  ,1  Jl  + 1 H-CIl*  (Ai.3  ( I,  J.l  , XFM.  JP,  NJ,  NF  ,2,21  ) 
T<Il,K2+lJi>=T(Il,K3+IJi)+Clt¥ALB(I,J,C,FM,JF,NJ,NF,2,2) 

S( I li I Jifl ) = S(I1  , I Jl  + 1 )+ClC*  ALd( I ,J,2 , FMP, JP,NJ,NF, 2,2) 

63  IJ2  = IAES ( I - J ) 

IF ( I J2 .GT . KFCUR)  GC  TO  5 

T(Il,IJ2+l)=T(Il,lJ2+i>+Clr»ALU(I,J,4,XFM,jF,N~,NF,2,2) 
T<I1,K3+1J2)=T<I1,K3+IJ2)  +C1C*AlB(I,J,5,FM,JP,NJ,NF  ,2,2) 
5(Il,IJ2+l>=S(Ii,IJ2+i)+ClC*ALU<I,J,6,FMP,jP,NJ,NF,2,2> 

5 CONTINUE 
3 CONTINUE 

CCMPATI EIL  IT Y PCUATIJNS  FOk  1=1 ■ 2 ,,,,,, 2KFOUR 
El=Cir/2. 

1 2 = K j.  + 1 
13=12+ K2-1 
CO  7 11=12,13 
I=Ii-Kl 
IS= I * * 2 

R ( T I » I1+K4-K1) -Oil 
T ( 1 1 , 1 1 +<4-K  D = - IS*C11 
T (II,  ID  =IS¥¥2*Clc 
IF (I.GT.KF  UUR)  GC  TO  82 
P (11,  I1+K3-K1)  = Q 1 1 
T(11,I1+K3-K1)=-IS¥C3+1./kP 
T (II,  I1-<1+1)=IS¥¥2*C8 
82  CO  9 J1=1,K1 
J=J1-1 

T(Il,v.l  + K2-l>  = T<Ii,Jl  + K3-i>-CK'MALB{I,J,l»WZ,JH,NJ,NW,l,l)  + 

1ALU (I , „,4 , WZ,JP , NJ  ,NW,  1,1 ) ) 

T (Il,Jl)=T  (Ii»Jl)-GiC*  (ALB(I»J»2  ,hZPP*jP*Nj,NW,i,l)  + 

1 ALB ( I , J , 5 ,WZPP , JP  ,NJ, NW, * ,1) ) 

SI  1 1,  Jl)  = S<Ii,Jl)-UluMAUU<I,J,3  , WZP, JP ,UJ ,NW, 1 , 1 ) + 

1 ALU ( I , J , 6 ,WZP, J F, NJ ,NW  ,1 , i ) ) 

IF (LN.EO.l)  GO  TC  9 

T(I^,J:  + K3-1)  = T(I1,J1  + K3-1)-E1*(AI?(I,J,1,NM,JF,UJ,N.I,1,D  + 

1 AL3(I,j,4,KP,JP,NJ,NW,1,1)) 

T < 1 1,  ul)  =T  (II,  Jl  ) - El*  ( ALo  ( I , J,2,  ETM,  JP,NJ  ,(JW,1 , 1)  ♦ 

1 ALB(I ,J»5»ETM,JF,NJ»NK,1,1) I 
S ( Il,Ji) =S (II, Jl )-ElMALO( I , J , 3 , WPF,JP,NJ»NW*l*i)  + 

1 A L ri  ( I , J , 6 » t*PP  , J F , NJ  » NW  » 1 , 1 ) ) 

G(  I1,NkH3)  =G  IL *NRPS)-El*  ( ( A L B ( I , J,1,WF,JF,NJ,NW,1,  1) 

1 + ALB ( I , J , 4 , rtf', JP  , N J , 

2 N W , 1,  1 ) ) * E T N ( J P , J 1 ) + (ALOt  1 ,J  ,<.*ETP»JP»NJ»NW,1,  1)  + ALB<i.,J,r.  ,E7P, 

3JP , NJ , NW , 1,1)) *WP(JP»Ji>+  (ALE (I,o,7,WnF,JP,NJ,NW,l, 1) +ALB(I,J,b, 

4 Is  M P , J F , N J , N A , 1 , 1 ) ) * WM  P ( JP  , J 1 ) ) 

I = T ♦ J 

tF  ( I J * .GT  .KFCUF ) GC  TO  83 

T (Ii,  Iw.I  + 1 ) = T ( II  *IJ1+1)"E1¥  ( A L L;  ( I ,J,l,cTM,JP,NJ,NW,  2,1)  ) 


/ & 

T(Il,Kc  + Idll=T(Il,K3  + IJl)-El*(ALaCI,Jt2,WM,JH,NJ,NHft,in 
S(Ii,lJl  + l)=S(Il,IJl  + l)-ElAfALb(I,J,3tW.vP»JP.NJ,NWtc,:» 

H3  IJ2  = 1 Af-S (I-J) 

IF ( IJZ .GT . KFCUR)  GO  TO  9 

T<I1,IJ2  + 1)=TCI1,Ij2+1>-E1*AI.G<I,J,4,ETM,JP,NJ,NW,2,i> 
T(Il,K3+IJ2)=T(Il»K3+IJ2)-Ei*ALB(I,J,5»WM,JP,NJ,NW,2,l) 

S ( II , 1*2*1  > =S(Il  ,1 J2-H  >-Ei*  Aid  <1  , J,6,  FiKP,  JP,NJ,NW,  2,1) 

9 CONTINUE 
7 CONTINUE 
RETURN 
END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

L, 

c 

<J 

c 

c 

c 

c 

c 


SUBPOLTINF  L,CUNDK(eS,BT,uG,IN,XNXX,lS,F'l,NJ,NW,NF,l_N,NfiFIS> 
COMNCN/FOUPIR/KF CUR , KF , K4 , K3 , K2 , K 1 
CCM;lCIWPRLSl/Wf'<2('C,5),ETM(2i)C,5>  ,WMH  (2.0,5) 

CCMNCN/PPESZ/WZ ( 2 c 0 , 5 ) , WZP  ( 2 C . , 5 ) , WZPP  ( cG  0 , 5 ) 

COMMON/PRES3/FNI  < 2?  ?,  8)  , XF  N ( 2 i ' , 6 ) ,FMP(2CP  , 8) 

COFdGN/GECf-V^R,  0C,P11 , Hi2,P22,Qll,Oi2,C22,Dll,D12,D22 
CCMMCN/FAC12/DL1  »UL2,UL3,CL4,[JA1  , 0A2  , D Ao,  0 A4 , Dd  2 , Or>  3 , DO  4,  XU  I , EX  XP 


COMMUN/FACTCR/C1 ,C2,C3,C4,C5,C6,C7,Ca,Cy,ClC,Cll,Cl2 
UIMENb  ION  OS  (Ml,  Ml  ) , BT  (Ml  ,K1)  ,BC»  < Pl.NRHSJ 
tCUNUAKY  C C NC IT  I 0 NS  8S*Z 5 + BF *Z=3 G 
L S=  1 FCR  SSI  W=NX=NXY=NX=N . 

k=KX=NXY=U=?. 

W = M X= V = NX  = 1 . 

W=^X-V=U=0. 

W=W*X=NXY=NX=P. 

W=W,X=NXY=U=T . 

W=W,X=V=NX=P. 

w =w , x= v = u=  c . 

NX=NXY=QX=PX=C  . 

NX Y = Q X=  W » X = l=  0 . 


L S = 2 FOR  SS2 
L S = 3 FCR  SS  3 
l S = 4 FCR  SS  4 
L S=  5 FOR  CC1 
l_S  = 6 FOR  CC  2 
L S = 7 F CP  CC2 
LS=  8 F CP  CCA 
LS  = 9 FOR  FREc.  EOGE 
LS=10  FOr  SYFf'ETRY 


LS  = 1-  FCR  AMISYFMETOY  NX  = P X = R = V= 0 . 


CALFA=(1.+XLAFC)*(1.+YlAFD)-XNI**2 

LL1=CD* ( 1.  + FhOX4CX  **2*XL AMU* (1 . +YL  AK0-XM**2) *12 ,/ ( Xh **2*0 ALFA)  ) 
CL2=LC*XM*  (l.+EX*EY*YLAF0*XLAKC*12  ./  < XH* * 2 * C AlF A ) ) 
CL3=EX*XLA.1CM1.  + YLAI''D)/CALFA 
DL4=-XrJl*EX*XLAMD/CALFA 
□Al=(l.t»LAMC)/(DALFA*EXXF) 

CA2=-XM/(DAlFA*EXXP) 

OA3=-(1.+YlAmC)*EX*XLAMD/OALFA 

DA4-XM*EY*YLflF'D/DALFA 

DE2=(1.*XLAMC)/(JALFA*EXXF) 

CE3=XNl*EX*XLAf'0/0ALFA 

CE4=- (1  . + XL A^C )*E Y *YLaMP/OAlFA 


K6  = fc  * K F OUR ♦ 2 


K4=4*KrCLRt2 


K2-3*KFClRt2 

K2=2*<FCIP 


K l=KFOUR  +1 

IF (LS .EQ. 11) LS=3 
CO  1 11=1, K6 
GG ( 11 , NRH  S ) =0. 

CO  1 -1=1 , KE 
ES  III  ,-i)  = ". 

E T ( 1 1 , J 1 ) = C . 

1 CON  TIME 

IF  (LS.&Q.  1'. ) GO  TO  888 
IF (LS.CT.8)  GO  TC  100 
DO  2 11=1,  Kl 

2 BT (11,11) =1. 

IF CLS.LS.4.CC.LS.GT.8)  GO  TO  1-? 
38®  J = 0 

00  3 I 1 = K 3 ,X  4 


o = J ♦ 1 

3 ns  ( n,  v.)  = l . 

10  ■>  IF  (Lo.tT,4)  Cu  TC  23? 
tT  (K*»K3)~i» 

K 3 1 = K ^ ♦ 1 


CO  4 I 1=  K 3 1 » K4 
OT (II, 11) =CL1 
t j T (I1,I1*K4-K3)=CL4 
IFtuS.EQ.l  .CF.LS.EC.3)  GU  TO  4 
1=  U-K3 

OT  (Il,Ii+Kl-KT)  = GT  (I1,I1-*n1-K2)-I+  + 2*C<;*DL3 

4 CONTTME 

2C0  IF ( LS . NE  . 1 r ) GO  TO  320 
£1=UL4/Dli+C10 
ES ( 1 , K 2) =CL1-DL4/011*011 
CO  5 v/i  = l i KFCUR 
JS=J1**2 

ES(l,Jl+i)=fiS(l,Jl+l)  + Ei*JS*VsZ(IN,Jl  + l> 
eT (1. - 1 + 1) =ET (1 , Jl+1) +E1* JS*hZP< I N, Ji+1) 

ET (1, K1+ J1 > =ET  <1 ,K1  + J1 ) +Clr  * JS#WZP( IN,  Jl  + 1) 

IF (LN.FQ.l)  GO  TC  5 

ES(l,2  1+i)=ES<l,JlU)+El*JS*fcM(IN,Jl  + i) 

ET (1, Jl+1 ) =ET(1 , Jl+1) +E1* JS+KPP< IN, Jl+1) 

EG(1.NFHS)=EGC1, NR HS)+E1*oS*M(IN, Jl+1) +WMP(IN, Jl+1) 

5 CON  TIME 

CO  12  11  = 2,  K1 
1 = 11-1 

OS  (II,  Il  + Kw-l)=t3S  (11,  Il+KE-1)  + UL1 
ES(I1»I1+K4-1)=BS  ( II , I1+K4-1  ) +DL 4 
CO  12  J 1=  1 , K 2 

BTCIl,Kl  + Jl)=nT(Il,Kl  + Jl)+ClC¥(ALE(I,Jl,2,WZP,IN,NJ»i'lW,l,l)  + 
1 ALB (I ,J1 ,5, fcZP , IN,NJ ,NW,1 ,1)  > 

12  CONTINUE 

3C  3 IF(LS.NE.S)  GO  TC  430 
E1  = LJL4/011  + C11 
ET (i,K2)=CLl-El 
dS (K3 , K3 > =Dll-£l 
E1=UL4/(U11*4R) 

E T ( 1 , : ) = - E 1 
LS(K3,l)=XNXX-wl 
OG<Kd.NRHS)=-XNX>*KZP«INt 1) 

E 1 = C L 4 / 13 1 1 * C 1 G 
CO  6 Jl=l, KFCUR 
JS= J1 **2 

ET (i,Jl+l)=ET(l,Jl+l) +E1+JS*WZ(IN,J1+1) 

ET  (Ki ,„1 +1) =LT (K2, Jl  + 1) +E1+JS+WZP  (IN , Jl+1) 

OS  ( < 7 ,Jl  + 1)  = OS  ( K 2, Jl  + 1)  + El‘rJS  + WZ  < IN, Jl+1) 

IF(LN.EQ.l)  GO  TC  6 

ET  (l,J  1+11=21(1, Jl+1) +E1*JS*WP(IN»J1  + 1) 
EG(i,NhHS)=Eo(l,NRHS)+El/c.*JS+MM(IN,Jl+l)**i 
OT  (K3,j1  + 1)=GT (K3,Jl+i)+£l*JS+WNP  (IN, Jl+1) 
dS(K?,Jl+l)=eS(K3,Jl+l)+El+JS*WM(IN,Ji+l) 

3G  ( K3  , NRH  3 ) = EG (KI.NRHS) +cl+JS  + WM ( IN, Jl  + 1) +WMP( IN, J*  +1) 
fc  CONTIME 

OC  *3  11=2, K1 
1=11-1 
IS= I*  *2 

bT ( II , II +K2-1) =DL1 
OT  ( II  , II  ) =-U*C9*0L2 
P T ( 1 1 ,11 +K4-1 ) = 0L4 
BS ( II 4K5-1 , I1+K4-1 ) =OL4 
BStI*«Kj-l,Il+K3-l)=LLl 


BS  < Ii  + KJ-1 , II)  =-IS*C9MDL2*2  .♦C0M1.-XNI)  ) + XNXX  ] ^ 

L+G  < I1  + K3  - 1 ,NKHS  ) =-XNXX*WZP  (IN,  Il  ) 

13  CONTINUE 
4Gr  I2-K141 

13-  Ic+KZ-l 

1 4 = K 3 - 1 

DC  7 11-12,12 
I = 1 1 - K 1 

15  - I*  * 2 

GO  TC  (21,22,23*  24. 21.26.27, 26,21. 3?) ,LS 

21  US (il ,K1 + 1 >=1. 

OT  ( II ♦ 14 , K1+  I) =1  . 

GO  TO  7 

22  liS  (II  ,K1  + 1 ) =1. 

os  di  + 14,  K4+i>=oe2 
00  8 21=1 , K 1 
J= Jl-1 

eS  < 11  + 14,  J1 ) =eS  t 11  + 14 , Jl) -C1C*  < ALO (I , J ,1 ,WZ,IN,NJ,NW,1, 1)  ♦ 
1AUB(I,J,4,WZ,IN,NJ,NW,1,1)) 

8 CONTINUE 

IF  (I.gT.KFOUR)  GC  TO  7 
CS(Ii+l4,K3+l)=BS(Ii+I4»K3+I)+0B3 
bS  <1^+14, I+l)=eS  (11+14,1+1 >-IS*C9+0d4+l./RR 
GO  TO  7 

23  BT  < 1 1 , Kl  + I ) = 1. 

OT  <11  + 14, K4  + I)=DE2 
IF (I.GT.KFCUR)  GC  TO  7 
BT (II +14, KZ+I) = D £ 3 
GO  TC  7 

24  ET ( 1 1 ,K1  + I ) =-IS* C9  + 3A2 
ET  ( 1 1 , K4+ I ) = CB2 

0S(I1+I4,K1+I) =-IS*C9* (DA2*2./((i.-XNl)*EXXP>) 
nS(Il+I4,K4+I)=D£2 
CO  9 J 1 = 1 , K 1 
J= Jl-1 

eS  <11  + 14, Jl) =dS  < 11+14, Jl)  - Cl"*  (AL  E( I, J ,1 ,WZ, IN,N J,NW, 1, 1)  ♦ 

1 ALB( I ,J ,4 ,WZ, IN ,NJ, NW,1, 1 > ) 

9 CONTINUE 

IF  U.LT.KFCUF)  GC  TO  7 
GT  (I1.K3  + I) = CB  3 
ES( 11+14, K2*I)=0E3 

PS  (11  + 14,  I +1)=BS  (11  + 14  ,1  + 1 ) -IS+U9*D34+*  ,/RR 

GO  TO  7 

cc  BS  < Il  ,K1  ♦ I ) =1. 

DS<I1+I-+.K4+I)=DE2 
IF  (l.GT.KFOUR)  GC  TO  7 
BS  (11+14, K2  + I»=0e3 
GO  TC  7 

27  CT < Il ,K1+I )=1. 

ET ( Il ♦ 14, K4  + I) =D  E2 
IF  ( I. G J . KFuUR)  GC  TO  7 
ET ( 11+ 14, K3* I) =0E3 
GO  TC  7 

28  OT(Il,kl  + I) = -IS*  C9*DA2 
dT ( 11 ,K4  + I ) =DB2 

BS(I.»+I4,K1  + I)=-  I3+C9*  (Dtfc+2,/((l  , 

BS( 11+14, K4+I)=DE2 


-XNI ) *EXXP) ) 


« 


IF  ( I.G  T.KFCUR)  GC  TO  7 
faT (Ii,K3*I)=CB3 
LJS  ( II  ♦ In  , K 3+  1)  =DP3 
GO  TO  7 

3:  GS  ( 1 1 » K1  ♦ I ) = 1 • 

BS  (I1*I4,KU  + I) = 0 G2 
DO  14  Jl=l * KFOUR 

BT  (11  4 14,  J 1 ♦ 1 ) = QT(  11*14,  JlU)-CirMALB(  I,  Jl,  2,  WZP,  I N,NJ,NW,  1,1)  + 
1ALU (I , J1 ♦ 5,WZP, I N,NJ,NW.1.1) > 

14  CON  TIME 

eS(Ii+I4,K3+I)=8S(Il+I4,K5+I)+033 
7 CONTINUE 
RETURN 
END 


I 


SUE  ROUTINE  INVERT(M»A,C.h.Nrtl,NH2,dcT»IXP,IOc.TJ 
DIMENSION  A (NFI 1 * NH 1) , C (N M2 ) , M (NM 2 ) 

CtT=i. 

IXF  = 0 
NN  = NA 

IF (NN. NE. II  GO  TC  803 
CET-A(ltl) 

A< 1 ,1) =1 ./A  (1,1  I 
GO  TO  ;nA 
30  3 00  ■)  C 1=1  , NN 
9 0 F ( I ) = - I 

DO  lAf  11=1, NN 
D=3 .DC 

DO  112  K=  1 » NN 
IF  ( FI  ( K ) ) lr',,lnU  ,112 
100  DO  11.  L = 1 , N N 

IF  < V < L > > 1 13,103, 110 

10  3 IF (AOS  (U ) - A E S ( A ( K , U ) ) ) 105,105,11  0 
1C  5 LO=L 

KD  = K 
0= A (K, L> 

BIG A=D 
11C  CONTINUE 
112  CONTINUE 

IF (D. tC. C . L J ) GO  TC  170 
GO  TO  138 
17G  WFITl(G,5'12) 

STOP 

3C  2 FORMAT  (/ , 5 X , "DE T ERMIN ANT=(  "/) 

18  ? N E V p = - F (LO) 

F ( L □ > = N (KC) 

F (KC)  = N E MP 
DO  11 A 1 = 1, NN 
C (1) = A (I, LC) 

A ( I , LD)  = A ( I, KD) 

11  A A ( I ,KC)  = C .CO 

A(KD,KL)=1.0F 
CO  115  J = 1 , NN 
115  A(KJ»„)=A(KD,J)/G 
LO  135  1 = 1, NN 
IF(I.EG.KD)  GO  TC  135 
CO  1 3 A J= 1 , NN 
TEMF=C(I)*A(KO,J) 

18  A A ( I , J) =A ( I , JJ-TEFP 
185  CONTINUE 

IF (IDET.NE.l)  GO  TC  1 A 0 
GET = GET* 9 IGA 
IF  ( Kl)  .NE  . LC)  CE  T = -OET 
:>29  IF  ( AeS  (DE  T ) . IT  . 1 .E+ir)  GO  TO  to0 
LLT=C-T/1.E+10 

IXF=lXF+in 
GO  TC  F29 

•33  0 IF  (AdS  (OET  ) .GT  . 1 ) GO  TO  lAn 

C c.  T = 0 r T*1  . £ ♦ 13 
IXr=IXF-10 
1AC  CONTINUE 

DO  2 1 = 1, NN 


'I 


1=0 

15  p l - L ♦ 1 

IF  < M ( L )-I  I 150 »1 

16  C Kl)=HIII 

F ( I ) = I 

DO  20 n J=1»NN 

TEMP=P  (L.  J> 

A (L, J)  = A (ltJ» 
2CC  A(I,J)=TEf'F 

30  a return 

END 


t 


1H 

1 ,15p 


A 


SUCROIUNE  YNY<N1,A,U,C,Nc,L1,l2,L3,T> 
DIMENSION  A ( LI  ♦ L 2 ) »B(Ll,Li),C(Li»L2)  , T ( L 3 ) 
IF ( N2 . £Q • 1 ) GO  TC  IOC 
DO  11  1=1, N1 
CO  U J=1,N2 
T £ MP  = n . 

DO  2u  K=1,N1 

20  TEilP  = Tti'iP  + 9(I,K)  *C  ( K , J ) 

10  T ( J ) = T E Mp 

DO  30  o=l, N2 
30  A ( I , J) =T ( J) 

11  CONTIME 
RETURN 

IOC  DO  111  1 = 1, M 
T EM P=  C . 

CO  120  K = 1 * N1 

120  TEMP=TEMP+E<I,K) *C<K,1) 

111  T ( I ) = T EMP 

00  13l  1 = 1,  M 
133  A<I,11=T(I» 

RETURN 
E 


SUJ  ROUTINE  Y S Y N Y (N2,Nl,A,d,C»D,N3,Ll,L2,L3,LA,T  ) 
CIMbNSlON  A (LI i L 3)  , B (LI , L3>  , C (LI  , L2)  , U (L2 , L3>  , T (LA) 
IF  (Mi  . *IQ  • 1 ) GO  TO  ICC 
00  11  1 = 1, N1 
CO  If.  v=  1 , N 3 
TEMP=C . 

DO  2 C K=1,N2 

2 ° TEMP=TEMP+G(I,K)  *0  (K,J) 
lH  T (J)=l-  (I,J)-TEMP 
00  3C  J=i,N3 
3 C A ( I » J J =T (J) 

11  CONTINUE 
RETURN 

1j  u 00  111  1 = 1, N1 

T c MP= C . 

00  120  K = 1 , N 2 

12  C TEMP=TcMP  + C (I,K)  *0  IK, 1) 

111  T (I»=0 (I , 1 ) -TEMP 

CO  1 2 C 1=1,  N1 

13  0 A (1 ,1)  =T  ( I) 

RETURN 

END 


SUBROUTINE  PCTEKS(IO:.T,N«HS,N.AXN,AP,DP,CP,GP,PR,XP,C,MT,Tl, 
1 V1.NAX2, IXFM,CETM,XNXX,LN,NJ»NW , NF ) 

M AX2=MAXN*MAXN 

COMMON /CINTG/NEQF0T,MI(5CC> 

COMMON /CO  I SK/ 121  < 3 1 1 » , 122  < 52  1 ) 

Di MENS  ION  AF (MAX  N , MAXN  ) ,RF (M^XN,  MAXN) ,CP ( MAXN.MAXN) 
DIMENSION  PR (MAX  N ,MAXN) ,GP (MAXN, URHSI  ,XP(MAXN,NPHS) 
DIMENSION  T 1 (MAX  N)  ,C (MAXN)  » M T ( M A X N ) * V 1 (MAX2) 

cQuivalence  (af  (1,1)  ,vi(in 

ixfm=c 

CE  TM  = 1 . 

DO  IOC  1=1. NEOPOT 

CALL  A ICG ( I.MAXN ,CP, BP , AP.GP .NRHS.XNXX , LN, NJ.NN, NF) 

N = MI  ( I ) 

iF(I.EC.l)  GG  TO  808 
NMIN1  = K(I-1> 

308  IF (I . tC.NEQFCT)  GO  TO  E99 
NPLtlal =M1 (1+1) 

99°  CONTINLE 

IF(I.cC.l)  GO  TO  12 

CALL  VSYMYtNMINl.N.eP.EP, CP, PR,N  ♦ MAXN, MAXN, MAXN, MAXN, Tl) 

12  CALL  INVERT (N, OP, C,MT, MAXN, MAXN, UET,IXF,IO^T) 

IF(ICET.NE.l)  GO  TC  6 4 C 
CETM=LLT*DE TM 
IXPM= IXP+IXFM 

IF (ABc (DLTM) .LT , 1.E+1C ) GO  TC  630 
DETM=CETM/1 . t ♦ 1 0 
I X F M = I X P M ♦ 1 T 
GO  TC  c 4 0 

63"  IF < A d S ( D E T M ) . G T . 1 . E - 1 r ) GO  TC  64" 

DETM=CETM*1.E+10 

ixpm=i>pm-io 

oh  Ci  CONTINLE 

IF(I.FC.NECFCT)  GO  TC  1C2 

CALL  VMY (N, FR, BP , AF , NPLUSl , MAXN, MAXN, MAXN, Tl) 

CALL  /FRITS  (21,PR,N',NPLUol,MAXN,MAXN,I,MAX2,Vll 
i: 2 IF  ( I . EC, 1 ) GC  TO  32 

CA  <_L  YSYMY  (NM.IN1  ,N,XP,  GP,  CP,  XP.NRHS,  MAXN,  MAXN,  NRHS,  MAXN  ,T1) 
CALL  >MY(N,XF,BP,XF,NPHS,MAXN,NRHS,MAXN,T1) 

GO  TO  42 

3?  CAlL  Yf Y (N , XP.OP ,GF,NRhS, MAXN, NRHS, MAXN  , Tl ) 

42  CALL  XFRITE  (22,  XF,N,  NPHS,  NAXN.NRHS,  I ,M,AXN,  Tl) 

U3  CONTINLE 

MEJPCT  = NECFCT-1 
DC  23.  K=l,MEQPOT 
NK=NEuFOT-K 
NM INI =MI  INK) 

N-MI(NKfi) 

CALL  XRtAC(2l,FK,NMlNl,N,NAXN,NAXN,NK,MAX?,Vl) 

CAlL  X R ” A C (22, GP  ,NM1 N1 ,NRHS, MAXN , NRHS, NK, MAXN, T 1 ) 

CALL  YSYMY<N,NMIM,XP,CP,FR,XF,NRFS,MAXN»MAXN,NRHS,MAXN,T1) 
CALL  XNRI TE (22 , XF, NMINI.NKHS, MAXN ,NRHS ,NK, MAXN ,T1  ) 

23°  CONTINLE 
RETURN 
END 


SLcKullINE  XKcAD (NO, A ,L1,l2»K1,M2 . INDt M3, VV) 

COVMCN/CO  IS</I21  (5  Cl)  , 122 (501) 

DIMENSION  a (Ml, M2)  ,VV  (M3) 

FECORO  IND  CF  uIRECT  ACLESo  DATA  SET  NO  IS  R;_AD  AND  ALLOCATED 
BY  FGRS  1MC  l 1 * L 2 PORTION  OF  MATRIX  A 
L3=L1*L2 

CALL  KEAuMS(NO»M*L3,  IND) 

KL  = C 

DO  1 NR0K=1»L1 

CO  U NC0L=1»L2 
Kl  = KL  ♦ 1 

A (NROW  ,NCCL) =VV C KL ) 
in  CON'TINLE 
RETURN 
END 


SUBROUTINE.  X*RITE(N0tA»LltL2*F,I«M2*IND»M'?tVV> 
CQFiTCN/CUlSK/IC*  ( 5 ^ i.  * » 122  (503) 

DIMENSION  A<Fl,Mcl  .W<F3)  Arrrce 

C L1*L2  FCft T’lON  OF  FA  T R 1 X A IS  WRITTEN  BY  ROWS  ON  DIRECT  ACC  EoS 

C CATA  SET  ND  IN  RECCFC  INC 
KL  = 0 

DO  10  NROW=l»Ll 
CO  1C  NCOL=l»L2 
KL=KL+1 

\]\l (KL)=A(NRCW,NCCL) 

10  CONTINUE 

CALL  WRITMS(ND»VViKLiIN0i-1) 

RETURN 

END 


SUBROUTINE  =uTSN  <F  CT  ,S  TRY,  STM,  IP,  1 SY,  ISA.  XNXX)  r>  Q 

C F CT  - FOTFNTIAL  ENtRGY  y 

C STkY  - UNIT  £NC  SHORTENING  FCR  Y = U . 

U STRA  - AVERAGE  UMT  CNC  SHORTENING 

U I F= 1 FCR  CALCULATE  F 0 T 

C IS  Y = 1 FCR  CALCULATE  c T R Y 

C I S A= 1 FCR  CALCULATE  STRA 

COMMON /FOUR  lk/KF  CUF,Kfc,K4,K3,K2,Kl 

CGMMCN/GFCP/RR,OC,P11.F12tH22,C11 ,0i2»CZ2,Ull,ul2,D22 
CGPPON/PRES1/WP  ( 2 C n » 5 ) ,cTM(2uf  ,5  ) ,W.1P(2on,5) 

COMMON/P RES2/W2 (20C,5) ^WZPtZCC,^*  ,WZPP  (20C ,5) 

CCMKCN/PRES3/FN ( 2: r ,fl)  tXFM (2CC » 8)  »FMP (2CC  » 8) 

COMMON /F  ACTCS/C1  ,C2,C3,C4,C5,Co,  C7,Ca,C'3,CiG,ul^,C12 
CCMPCN/FACT  2/DL  1 *DL2*DL3»DL4,CAl  »CA2*DA3,0AA,UE2*0u3»D3A*XM»EXXP 
COPMON/F  AC  T3/CL5  tXL*  XH 
COPMUN/CINTG/NEUFOT , M I (5b f) 

CO  M PC  N /FICFR/DEL  TA  » A L 1 * G 4 1 »AL2»B  T 2,GA2 

CONHON'/XXLOAO/XPRES 

CEl=ClC/2. 

CE2-CE**2/2. 

CE  3 = C E / ( (l.-XNI) *EXXFJ 
CE  4 = C° *DD* (l.-XNI) 

FO  T = 0 . 

STR Y-o  . 

S T R A = 0 . 

DO  1C  I1=1,NECPCT 
£7  = 1. 

IFdl.tQ.l.OR.Il  . tG.iMECPOT)  E7  = C . 5 
E1=-C11*ETM(I1,1)-1./RR*WM(I1,1> 

E 2 = 0 . 

hO  li  jl=l,KFOUR 
JS=J1 * *2 

E2=JS*kM(Il,Jl+l»MWM(Il,JH-l)+2.+RZ<Il,Jl  + l>MF.2 
11  CONTINUE 

E1=E1+CE1*E2 
IF (IF.NE.l)  GO  TC  lCu 

FEi  = bB2/DU**2*£l*,’2HJLmTf1<ll,l>**c-XNXX*WMP(Il,l)*(khP(il.l)  + 
12. *W7F  (II » 1 ) ) 

FE1  = P=.1-2.*XFRES"WP(I1,1) 

1GC  El=El*CA2/Ci;*DA 2*ETM  ( 11,1) 

IF(IGY.NE.i)  GO  TO  llr 
P S Y = E 1 

11C  IF  ( ISA  .NE  . 1)  GC  TO  12  8 

PSA=>_1-WHF  ( lit  1)  * ( WMP  (Il»l)+2.*W2F(Il,l)>/2. 

120  IF ( ISY  .NE  .1)  GC  TO  13r 
El  = 1 . 

DO  lb  o 1 = 1 * K 1 
CO  18  J2=l » K 1 

El  = WMF  (1 1 , J1 1*  (WPP  (1 1, J2>  *2.*  NZP ( II ♦ J2) ) + E1 
18  CONTINUE 

FS Y = PS K-El/2  . 
il  = 0. 

E 2 = 0 . 

00  u Jl=  1 » Kf-OLR 
JS  = J 1 * * 2 

lI-.TN  (I1,J1*1) +61 
E2=c.2*„S*N,JCIlf  Jl  + 1) 


12  CON T IN LE  2>  A 

PSY=PSY+nA3*tl-Dfl4*C9*E2 

E1=0. 

E2  = C . 

DO  13  Jl  = l *K  2 
wS=Jl«*2 

El=El+XF9  (II  ,J1» 

E2=E2+cS*FH(Ilf Jl) 

13  CON  r I ME 

FSY=PSY+0A2*E1-0A1*C9*E2 

STRY=STPY*FSY*F7 

13  G IFdSfl.NE.il  GO  TO  140 

E1  = 0. 

00  14  J1=1,KFCUP 

14  E 1 = ti  + K1F ( 11 * Ji  + 1) * (WMF (II , J 1 +1 ) ♦ 2 . * WZF (1 1 * Ja+ 1 > I 
PSA=PSA-Fl/4  . 

STkA=ilRA+FEfl*E7 

14  0 IF  (IP.NE.l)  GO  TC.  10 

E 1 = 3 . 

62  = 9. 

E 3 = C • 

E 4 = 0 . 

E 5 = r . 

DO  15  Jl=l,KFOUR 
JS  = Jl " *2 
v,S2  = JS”  2 

El  = Ei  + MPdl,JH-li'<WMF(ritJi*lM2.*WZF(Il,Jl*l)> 

E2  = E2  + .S2*NM(II,Jim**2 

E3  = E2+JS*MF<Il,Jin>**2 

E4  = E4  *ETM ( II ,J1  + 1) **2 

E5  = E-3*  JS*  WF.  (II  , Jl*l)  *ETM(I1,  Jl+1  ) 

15  CCNTINLE 

PE1-FE  1-*NXX*E1/  2.  ♦ CM*0LS*E2+CE  4*E3  ♦0Ll4‘£4/2  . -C9*CL2*£5 
E 1 = n . 

E 2 = 0 . 

c.3  = p . 

E 4=  C • 

UO  16  ol  = 1 » K 2 
JS=  J 1 ■*  * 2 
.32= JS**2 

El= JS2*FM  ( II  .Jl  > **2  + cl 
E 2 = E2  ♦ JS  * F F'F (II, Jl)»*2 
E3  = E3*>FM<I1»J1) **  2 
E4=£4+.S*FP(11,J1)*XFM(I1,J1) 

16  CONTIME 

PEl  = PEl*-CE2»CAl*El+CE3*E2*D3  2,'FJ/2.-0A2*C9*E4 
F0T=FC  T+F£1*F7 
10  CCNTINLE 

IF(iF.NE.l)  GO  TO  1 5 C 
F0T  = rCT*3.14l59*RR,,DFLTA 
150  IF(ISY.NE.l)  GO  TO  lbr 
LTr’Y  = Cfll*XNXX-ST  Sy/XL 
160  IF  (ISt  .t)E  .1)  Go  TO  17C 
STR A = C A 1 * X N X X - S T RA/XL 
173  CCNTINLE 
F E T U P N 
END 


SU  BROL 1 I NE  CCE.F  F <FX,EY  , XL  MID,  YL A NO , RHU X , RFIOY , E l AS  > 

CCMMCN/Gc.GM/RR,OC*Hli,H12,Fl22*Qll»Ql2, 022*011, 012, 022 

CCMI1CN/FACT2/OL1  ,JL2,DL3,l  L4,LA1  , C A2  , l)  A 3 , D A 4 , 0 6 2 , 0U3  , OE  * , XN  I , E X XP 

COMiION/CINTG/NlQFOT.MI  (5CD 

CCMMCR/FICM/DFLTA,ALi,GAi,AL2,3T2,GA2 

COKMCN/FOLRIR/KFCUF.Kb ,K4,K3,K2,  K1 

COMMON/F  ACT3/CL5  » X L , XH 

K6=G*KF0UR+2 

K4=4*KF0UR+2 

K3  = 3*KF0UR+2 

K2=2*KF0LR 

Kl-KFCIR+l 

XN2  = XM**2 

Xh2=XF»*2 

CALFA = (1  . + XL AM Cl  * ( 1 . + YlAMG) -XN2 
0D=ELAS*XH**3/(i2.*(i.-XN2) ) 

EXXP=tLAS*XH/(l.-XN2» 

HI 1 = 1 . 4RHCX  + 12. * EX**?*  XLAKU* (1.  + YLAMD-XN2) / (XH 2*0  ALFA) 

H22=l . «RHCY*l2. * EY**d*YlANO*  ( 1 . ♦ X L A.1D-XN2  > / < XH  2 * 0 A..  FA  > 
H12=i.*l2.*XNI*EX*EY*XlAML*YLAMD/(Xri2*CALFA) 

C11  = XM*EX*XLAM0/DALFA 
G22  = XM*EY*YLANU/0ALFA 

G12  = -:.3,‘«(l.*YLAMD)*EX*XLAMC+a.«-XLAML)*EY*YLAMU)/0MLFA 
C11=(1.*XLAMC)/(CALFA*EXXP> 

022= ( 1.* YLAKC) / ( CALFA *E XX P) 

012=1(1.  + XLA  MO  )*(1.*YLAMD)-XM)/(CAl  FA*  lX  X P*  ( 1 . - XN 1 ) ) 

LL1=CC*H11 

CL2  = Db*XNI  Ml.  + < EX*E  Y*XlANC*YLAMU*12  . ) / (XH24DALFA>  ) 

DL?  = EX*XLAMDM1.  iYLAMU)/DALFA 
CLA  = -XM*EX*  XL  AML/C  ALFA 
OL5=LL*H22 

DA1=  (1  .♦  YL  Ar1C»  / ( CALFA*  EXXF) 

0A2=-XM/  (LALFA*EXXP) 

CA3  = -C.+YLAK0)  *EX*XLAMC/LALFA 
0 A •,  = X N Y * YLAMO/OwlF  a 

OB  2=  ( 1 . ♦ XL  AND  / ( CALFA*EXXF) 

C32  = XM*EX*XLAM  0/0  AlFA 
CBA  = -(1.  + XLAFC)  *EY*YLANC/LAlFA 
MI (1) = 2*KE 
MI  (MGPOT  ) = 2*  K 6 
NEG1=NEQPCT-1 
CO  1C  11=2, NEQ1 
MI ( II > =K6 
1?  COMIME 

DELTA=XL/  (NtCPCT-1) 

AL1=-1./(2.*CELTA) 

GAl=l./(2.*0cLTA) 

AL  2 =1 ./(0ELIA**2  ) 

£T?  = -2  ./DELT  A**2 


5UCKCL11NE  CCEFNN(KNN) 

COdHUN/GEOM/^Pf  OD,  HU  ,H2,H22,U11,Q12,Q2Z,l»11,l)12,[;22 
CcPMCN/FACTCR/C  1.C2,  C..:  ,C<. , CEtUb,C7,C6,L9,Clr,C^l.C12 
C.9  = (NNN/RK  ) **Z 
ci  = uu*m 
C2=2.*C0*H12*C9 
Co=2. *Clt*C9 
C*  = L)D»F22*C9**2 
C5-m:/C)ll*C9 
C6  = 1 . / (RR*C1 1)  *C9 
C7  =09**2/  ( 2 • *D1 1 ) 

Cfi  = C2d*i,9**2 
ClC=C9/2. 

L11=2.*C9*D12 

C12=D22*C9**2 

RETURN 

END 
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